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A numerical  model  for  the  prediction  of  nearshore  hydrodynamics  of  an  inlet-beach 
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ing perturbation  technique  is  applied  to  arrive  at  the  same  wave  equation.  The  application 
of  3rd-order  wave  equation  is  restricted  to  small  incident  wave  angle  and  slowly  varying 
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CHAPTER  1 
INTRODUCTION 


1.1  Statement  of  Problem 

Tidal  and  river  inlets  are  prominent  coastal  features.  They  play  a vital  role  in  the 
evolution  process  of  coastlines.  Along  the  east  coast  and  gulf  coast  of  the  United  States, 
inlets  are  numerous  because  of  the  extensive  barrier  island  system.  In  these  regions,  the 
rules  of  inlets  are  particularly  significant  owing  to  the  heavy  development  pressure  and 
the  competing  interests  of  various  human  activities.  To  maintain  navigation  and  water 
commerce,  for  instance,  often  requires  inlets  to  be  stabilized  with  man-made  structures. 
The  disruption  of  littoral  processes  owing  to  the  presence  of  inlet  structures,  on  the  other 
hand,  poses  a major  threat  to  down  drift  beaches,  thus  diminishing  the  value  of  recreational 
beaches  and  endangering  waterfront  structures.  Many  coastal  communities  spent  millions 
of  dollars  installing  coastal  protective  structures  or  simply  transporting  sand  from  one  place 
to  another.  However,  these  measures  do  not  always  render  satisfactory  results  owing  partly 
to  our  poor  understanding  on  complicated  interactions  of  an  inlet-beach  system. 

The  hydrographic  and  shoreline  changes  of  an  inlet-beach  system  are  the  consequence 
of  sediment  movement  which,  in  turn,  is  caused  by  the  hydrodynamic  forces  exerted  on  the 
system.  Therefore,  to  better  predict  the  behavior  of  inlet-beach  interaction  one  must  first 
be  able  to  predict  the  hydrodynamics. 

Physically,  the  problem  in  hand  is  a complicated  one.  Both  wave  and  current  field  are 
modified  by  the  nearshore  topography  and  the  inlet  geometry.  In  addition,  the  spatially 
non-uniform  current  will  interact  with  the  unsteady  flow  field  induced  by  waves. 

Mathematically,  there  is  no  unified  wave  theory  in  the  first  place  capable  of  describing 
waves  traveling  through  large  areas  of  varied  depth.  Furthermore,  the  current  field  that  we 
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are  dealing  with  has  a different  scale  of  variation  from  that  of  the  local  wave  field.  The 
nonlinear  interaction  between  the  two  systems  is  expected  to  be  important  as  well  as  the 
non-linear  behavior  of  the  wave  field. 

1.2  Peist  Study 

In  the  absence  of  an  inlet,  the  water  motion  in  the  nearshore  zone  is  primarily  induced 
by  waves  from  offshore  as  modified  by  the  bottom  topography  and/or  shore  structure. 
The  influence  of  bottom  topography  on  wave  direction,  height  and  length  has  long  been 
recognized  and  as  early  as  in  the  1940s  analytical  methods  were  already  being  developed 
to  solved  the  so  called  wave  refraction  problem  (Johnson  et  al.  1948;  Arther  et  al.  1952) 
by  tracing  wave  rays  based  on  the  principle  of  geometrical  optics.  Since  then,  significant 
progress  has  been  made  in  treating  the  shallow  water  wave  transformation  phenomenon. 
Arthur  (1950)  included  the  effect  of  current  into  the  ray  method.  Dobson  (1967)  developed 
one  of  the  first  computer  codes  to  trace  wave  ray  over  an  irregular  bottom.  Skovgaard  et 
al.  (1975)  extended  the  ray  method  to  include  the  effects  of  wave  energy  dissipation  due 
to  bottom  friction.  Noda  et  al.  (1974)  developed  a finite  difference  scheme  based  on  the 
conservations  of  wave  number  and  wave  energy,  that  compute  spatial  distributions  of  wave 
number  and  wave  energy  in  shallow  water.  His  method  was  extended  by  Wang  and  Yang 
(1981)  and  Chen  and  Wang  (1983)  to  wave  spectral  transformation  and  the  inclusion  of 
non-stationarity,  thus  enabling  the  computation  of  spatial  and  temporal  distributions  of 
random  waves.  Birkemier  and  Dalrymple  (1976)  and  Ebersole  and  Dalrymple  (1979),  on 
the  other  hand,  expanded  Noda’s  method  by  including  the  depth  averaged  continuity  and 
momentum  equations  to  computer  wave  induced  nearshore  circulations.  In  the  presence  of 
an  inlet,  the  nearshore  wave  field  is  further  modified  due  to  the  interactions  with  a spatially 
non-uniform  current  and  the  diffraction  effect  of  inlet  jetty  structures. 

Traditionally,  the  problems  of  wave  refraction  and  diffraction  were  treated  separately. 
Penny  and  Price  (1944,1952)  were  among  the  first  to  apply  Sommerfeld’s  (1896)  solution 
to  wave  diffraction  around  a sharp  obstacle  such  as  a semi-infinite  thin  breakwater.  Other 
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solutions  of  similar  nature  were  obtained  for  different  geometries  ( Wiegel,1962;  Memos, 1976; 
Monfetusco,1968;  Goda  et  al.,1971). 

Only  recently,  Berkhoff  (1972)  developed  what  is  now  referred  to  as  the  mild  slope 
equation  that  governs  the  combined  refraction-  diffraction  phenomenon.  Various  solutions 
became  available  on  wave  field  modified  by  different  beach-structure  configurations  such  as 
waves  in  the  vicinity  of  a thin  groin  (Liu  and  Lozano, 1979;  Liu, 1982),  edge  waves  by  barred- 
beach  topography  (Kirby  et  ah, 1981),  and  waves  over  circular  shoal  (Radder,  1979),  etc. 
Booij  (1981)  further  introduced  the  ambient  current  element  into  the  basic  wave  transfor- 
mation equation  and  derived  the  governing  equation  on  variational  principle.  The  resulting 
equation  is  a hyperbolic  type  and  is,  in  general,  difficult  to  solve.  Booij ’s  equation  reduces 
to  a parabolic  type  similar  to  Berkhoff’s  if  the  mild  slope  assumption  is  invoked.  A number 
of  simple  examples  were  given  by  him.  Kirby  (1983)  systematically  studied  the  current- wave 
interaction  problem  and  proposed  a weekly  non-linear  model  appropriate  to  thr  Stokes  wave 
expansion. 

The  present  study  extends  the  Booij ’s  and  the  Kirby’s  model  to  the  problem  of  inter- 
action between  water  waves  and  inlet  flows. 

1.3  Scope  of  Study 

The  primary  purpose  of  this  study  is  to  develop  a numerical  model  for  the  prediction  of 
wave  characteristics  and  current  field  of  an  inlet  beach  system.  Specifically,  the  mathemati- 
cal formulation  of  the  problem  was  performed  first.  This  task  entailed  the  establishment  of 
the  governing  equations  sufficient  to  solve  the  variables  of  interest.  Out  of  necessity,  simpli- 
fied assumption  are  made  such  that  the  problem  is  amenable  to  the  present  computational 
capability. 

First  of  all,  the  fluid  is  treated  as  incompressible.  Second,  the  surface  shear  stress  eis 
well  as  part  of  the  turbulent  shear  stress  (normal  component)  is  also  neglected.  Third, 
also  one  of  the  most  major  ones  is  that  the  irrotationality  is  assumed  in  the  wave  equation 
allowing  the  existence  of  a wave  velocity  potential  and  a current  velocity  potential,  although 
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such  an  assumption  is  in  contradiction  to  both  physical  reality  and  mathematical  results 
of  the  combined  velocity  field.  The  equations  governing  the  mean  flow  including  both 
continuity  and  momentum  equations  are  depth  averaged.  The  equation  governing  the  wave 
motion  with  due  considerations  of  refraction,  diffraction  and  current  interaction  is  quasi- 
three  dimensional  owing  to  the  velocity  potentieil  assumption.  Two  parallel  methods  were 
used  to  derived  the  wave  energy  equation;  one  based  on  variational  principle  and  the  other 
on  perturbation  method.  By  introducing  proper  scaling  factors,  the  wave  energy  equation 
is  further  simplified  by  containing  terms  of  only  the  same  order  of  magnitude.  Other 
relationships  including  dispersion  relation,  radiation  stresses,  bottom  shears  and  lateral 
mixing  compatible  to  the  system  under  consideration  are  also  formulated. 

A numerical  scheme  is  developed  to  solve  the  governing  equations  by  double  sweep 
marching  method.  The  scheme  is  implicit  which  allows  for  a large  time  step,  especially  the 
two  time  level  fully  alternating  direction  implicit  scheme  utilized  by  Sheng  (1981). 

Attempts  to  verify  the  numerical  model  are  made  by  comparing  the  present  model 
with  known  analytical  and  numerical  solutions  of  simple  geometries  and  with  laboratory 
experiments  carried  out  in  the  Coastal  and  Oceanographic  Engineering  Laboratory.  Due  to 
wave  basin  limitations,  the  experiments  were  carried  out  only  for  plane  beach  with  an  inlet 
normal  to  the  contours. 


CHAPTER  2 

THEORETICAL  EQUATIONS  OF  THE  MODEL 


2.1  Introduction 

In  order  to  solve  the  flow  field  in  a specified  region,  the  governing  equations  with  proper 
initial  and  boundary  conditions  must  be  provided.  Furthermore  whether  a unique  solution 
exists  in  the  domain  of  interest  or  whether  mathematical  singuleirities  will  appear  that  will 
cause  the  equations  failing  to  give  sensible  solutions  should  also  be  examined. 

In  this  chapter,  the  theoretical  development  of  the  governing  equations  for  both  the 
mean  flow  field  and  the  wave  field  is  given  in  detail. 

There  are  five  basic  equations  employed  in  the  computation  including  four  flow  equa- 
tions: one  continuity  equation,  two  horizontal  momentum  equations  and  one  wave  equation, 
and  one  equation  relating  wave  number  and  wave  period. 

2.2  Governing  Equations 

The  basic  equations,  as  mentioned  before,  are  based  upon  the  principle  of  conservation 
of  mass,  conservation  of  momentum,  and  conservation  of  wave  energy  and  upon  the  gravita- 
tional wave  dispersion  characteristics.  The  wave  equation  under  consideration  includes  the 
effects  of  wave  refraction,  diffraction  and  shoaling  during  transformation  from  deep  water 
to  shallow  water  regions.  These  equations,  together  with  initial  and  boundary  conditions, 
give  a mathematical  representation  of  the  physical  problem  of  interest. 

Since  the  problem  of  interest  entails  waves  propagating  towards  shore  over  a complex 
bathymetry,  wave  refraction,  diffraction  and  breaking  are  expected  to  occur  as  well  as 
current-wave  interaction  due  to  inlet  discharge.  Figure  2.1  depicts  the  physical  system  being 
studied.  The  problem  has  two  apparent  time  scales:  a short  time  scale  associated  with  the 
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Inlet 
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Figure  2.1:  Outline  of  Studied  Area 
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oscillation  of  the  incident  waves,  and  a larger  time  scale  over  which  the  characteristics  of 
the  incident  wave  and  mean  flow,  such  as  wave  height,  incident  wave  angle  and  current  are 
perceptibly  modified.  The  larger  time  scale  may  also  include  the  effect  of  changes  in  wind, 
tides  and  topography.  Since  our  attention  here  is  towards  mean  quantities,  many  of  them 
may  be  averaged  over  a wave  period  to  remove  the  direct  effect  of  oscillation.  In  addition, 
since  we  are  mainly  interested  in  the  net  transport  quantities  the  knowledge  in  the  detciiled 
structure  of  the  velocity  distribution  over  depth  is  not  needed  and  the  equations  may  be 
averaged  over  depth,  reducing  the  entire  problem  to  two  dimensional  in  the  horizontal  plane. 

In  the  derivation  of  wave  energy  conservation  equation,  due  to  the  extreme  complexity 
in  mathematics,  some  restrictions  have  to  be  imposed  on  the  orders  of  magnitude  of  wave 
height,  bottom  slope  and  current  strength. 

In  the  model,  the  quantities  to  be  determined  are  as  follows:  the  horizontal  components 
of  the  mean  current  including  wave-induced  current  and  divergent  current  from  the  inlet; 
the  local  wave  ^lmplitude  under  the  influence  of  refraction,  diffraction  and  breaking  ; the 
local  wave  angle  eind  wave  number  and  the  mean  water  surface  level. 


2.2.1  Depth- Averaged  Equation  of  Continuity 


A Cartesian  coordinate  system  is  defined  with  x in  the  onshore  direction,  normal  to  the 
co8istal  line,  y in  the  longshore  direction  and  z verticeilly  upward,  as  shown  in  Fig. 2.1. 

The  summation  convention  adopted  here  states  that  in  Cartesian  coordinates  whenever 
the  same  letter  subscript  occurs  twice  in  a term,that  subscript  is  to  be  given  all  possible 
values  and  the  results  in  horizontal  space  add  together. 

The  continuity  equation  reads 


^ dpuj  ^ 

dt  dxi  dz 


(2.1) 


where  » = 1,2  represents  x,y  direction,  respectively,  p is  the  water  density  and  u,-  = {u,v} 
represents  the  {x,  t/}  components  of  velocity,respectively,  w is  the  velocity  component  in  the 


z direction. 
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In  application,  water  density  which  only  varies  slightly  is  considered  as  constant.  The 


continuity  equation  reduces  to 

dui  dw 
dxi  ^ dz 

Integrating  Eq.(2.2)  over  depth,  we  obtain 


(2.2) 


dw 


)dz  = 0 


(2.3) 


where  rj  is  free  surface  elevation,  and  h is  water  depth.  Using  Leibniz’s  rule  which  states 
that 


dx 


df  da  dd 

f{x,z)dz=  -^dz+ —f{x,a)  - -^f{x,l3) 

Jp{x)  j0(x)  dz  dx  ^ dx  ^ ' 


i-a(x)  Qf 
I0[x) 


Eq.(2.3)  can  be  written  as 

d /■’?  dr]  dh  , , , ^ 

^ I,  l-k  ^ I,  l-»=  0 (2.4) 

in  which  the  subscripts  after  vertical  stroke  are  at  which  the  value  being  calculated.  The  bot- 
tom boundary  condition  (BBC)  and  the  kinematic  free  surface  boundary  condition  (KFSBC) 
are  given  as 


u,—  + tu  = 0 (BBC) 


at  2 = —h{x,  y) 


(2.5) 


— = (KFSBC)  at  z = tj(i,  y,  t)  (2.6) 

Substituting  (2.5)  and  (2.6)  into  (2.4)  results  in 

dr)  d f')  , 

We  now  introduce 

«,•  = Ui  -h  u| 

v = v + n'  (2.8) 
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where  ?/,■  are  time  and  depth  independent  uniform  current  velocities,  ff  is  the  mean  surface 
elevation,  it  j and  rj'  are  the  fluctuating  components.  u\  can  be  separated  further  into  two 
parts: 

u'i  = uf  + «.•  (2.9) 

where  uf  and  u\  are  the  slowly  varying  wave  motion  and  the  fast  varying  turbulent  fluc- 
tuations, respectively.  The  properties  of  these  quantities,  after  taking  average  over  wave 
period,  are  expressed  as 

Ui  = Ui,  ^=ff 

^ = 0,  ^=0 

where  the  upper  bar  represents  the  time  average  over  one  wave  period 

- 1 

> = tL 

Therefore,  (2.7)  becomes,  after  taking  time  average 

dr}  d — ~d  7^ 

Defining 

. ^ I-hKdz  ^ !-h  <dz 

(h  + ff)  {h  + r}) 

the  final  form  of  the  mean  continuity  equation  is  arrived  at 

^ + (211) 

where  D = h + r)  is  total  depth  and  C/,- = 17, • is  total  mean  current  velocity  averaged 
over  one  wave  period  with  Ui  = {U,V}.  Equation  (2.11)  represents  conservation  of  mass 
per  unit  surface  area  under  the  assumption  of  water  density  being  constant  in  space  and 
time,  and  the  bottom  unchanged  with  time. 
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2.2.2  Depth- Averaged  Equations  of  Momentum 

The  horizontal  depth-averaged  momentum  equations  used  here  were  originally  derived 
by  Ebersole  and  Dalrymple(l979).  A brief  summary  was  given  here.  The  Navier-Stokes 
equation  in  the  horizontal  (x,y)  directions  takes  the  following  form: 


dui  dui  dui  \ dp  1 dTji  1 dr^i 


and  in  vertical  (z)  direction  it  is 


dw  dw  dw  1 dp  1 drj;i  1 dTzz 
dt  ^ dxj  dz  p dz  p dxj  p dz 

1 = 1,2,  y = 1,2 


(2.12) 


(2.13) 


where  p is  the  total  pressure  including  dynamic  and  static  pressures  and  r.-y  is  the  shear 
stress,  where  the  first  subscript  is  the  surface  that  the  stress  acts  on  and  the  second  subscript 
is  the  direction  of  shear  stress.  Introducing  the  continuity  equation  (2.2)  into  (2.12)  and 
(2.13)  we  have 

dui  duiUi  du.w  1 do  1 dr.-.-  1 dr^.- 

(2.14) 


^ dujUj  ^ dujW  _ 1 dp  1 dTji  1 dr^i 
dt  dxj  dz  p dxi  ^ p dxj  ^ p dz 


and 


dw  dwuj  dww  _ 1 dp  1 dr^j  1 dr, 

dt  dxj  dz 


+ 


pdz  ' p dxj  ' p dz  (21^) 

Integrating  over  depth  and  utilizing  the  Leibeniz’s  rule,  the  following  term  by  term  results 
are  obtained: 


i) 


n dui_dp  dr) 

J-H  dt 


duiUj  dr'? 


dr) 

da:,- 


dh 


r,^  - I-/.  ^ 


du{W 

iiij  y_^  ^ I”  I-/* 

Thus  the  left  hand  side  of  Eq.(2.14)  becomes 


LHS  = 


_d  n 
dt 


Lh  dxj  Lh  I"  ^ ~ I") 
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dh 


-«•  \-h  {Uj  l-fc  + W l_ft) 


Recognizing  the  last  two  brackets  of  LHS  are  the  KFSBC  and  BBC,  respectively,  we  have 


LHS 


a /•'J  , d n 


dz 


Again  introducing  the  definition  of 


u,-  = Ui  + uf  + u‘ 


into  the  above  equation  and  time  averaging  over  a wave  period  yields  term  by  term 

where  H7,-  and  D are  total  current  velocity  and  total  water  depth  respectively  as  defined 
earlier,  and 


“>  I-H  ° ^ + “S'* 

Noting  that  there  is  no  interaction  between  wave  orbital  velocity  and  turbulent  velocity, 
Phillips  (1977),  the  above  term  can  be  expanded  as  follows 


_d 

dx 


/n 

u\u\dz 
-h  * ^ 


_d 

dx 


~{U,UiD)  + ^ ufujdz  + 


— r u*uUz 
dxjj-h  ‘ ^ 


Introducing  Ui  = Ui  + ti,-  into  the  above  expression,  the  LHS  of  Eiq.(2.14)  becomes 
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5x,- 


(216) 

The  right  hand  side  of  (2.14)  is,  after  integrating  over  depth  and  averaging  over  time, 


1 r j 1 r dTji  ^ 1 n dr^i  , 

RHS  = - - / ^dz  + - / -^dz  + - / -^dz 
p J-k  OXi  p J-h  OXj  p J-h  oz 

here  water  density  is  considered  as  constant.  After  neglecting  the  second  term  which  is  the 
laterzJ  shear  stress  due  to  viscosity,  we  obtain 


1 ^ r j 1 I 1 , dh  1 , 1 

RHS  = - / pd2  + -p  ' + _p  _^  + _7-^.  _ 

p dxi  J-h  p OXi  P OXi  P P 


'^zi  — /i 


The  datum  of  pressure  is  chosen  eis  the  atmosphere  pressure,  so  that  at  free  surface,  p |,j 
is  zero.  The  mean  pressure  at  the  bottom  p \-h  consists  of  two  parts:  hydrostatic  pressure 
and  wave  induced  dynamic  pressure,  viz. 


p|-fc  = -P9{h  + r))+pd\-h 


= -pgD  + pd\-^ 


where  pd  is  the  dynamic  pressure.  Multiplying  the  above  equation  by  1 ^ we  have 

P 


1 I dh  15/  ^l2^  n 5»7  1 dh 


The  RHS  can  now  be  written  as 


1 ^ r j , 1 ^ 1 , dh  1 


-ya  (2.17) 

where 

^ si  = '''zi  |rj  is  time  averaged  surface  shear  stress 
and 

Tii  = T^i  |_fc  is  time  averaged  bottom  shear  stress. 
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Finally,  combining  (2.16)  with  (2.17),  the  time  averaged  and  depth  integrated  horizontal 
momentum  equation  is  given  by 


pdx,  ^ I-*  + /«  - 


or  rearranged  eis 


1 r»2c  - - r>l  , 1 I 1-  1-  d f'l  . . 

pgD  Oij  puiUjD]  H pd  |_;j  1 — r*,-  H — rj,-  — - — / u^u*-dz  (2.18) 

P oxi  p p dxj  J-h  ^ 

Collectively,  the  second  to  the  fifth  terms  on  the  RHS  are  known  as  the  radiation  stress  or 
the  excess  momentum  flux  due  to  the  presence  of  waves,  i.e.. 


^ ^ ^I-h  “ \p9^^^a  - pUiUjD  (2.19) 


where 


S-=l  1 ‘=^‘ 

1 0 

Considering  that  the  product  of  bottom  slope  and  the  mean  dynamic  pressure  at  the  bottom 
is  negligible,  we  obtain 


(2.20) 


after  assuming  ujuj-  is  independent  of  depth.  The  term  -pu\u*.  is  Reynolds  stresses  and 
their  normal  component  (when  t = j)  can  usually  be  neglected. 


Writing  in  component  form  yields 
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X-Direction 

+ §-^UVD)  = 


P dy 


Ddrt  1_  1_ 

7% +7“  “7"  (2-21) 


Y-Direction 


|(yD)+|-(t.yD)+|-(y^D)  = -,D|5- 


p dx 


P dy 


Ddrt  1_  1_ 

7 ^■^7''' ■7“'  (2-22) 

where  Yi  — -pu\v}-  is  the  lateral  friction  due  to  Reynolds  stress.  In  the  model,  r^,-  will  not 
be  considered.  The  possible  formations  of  rj,-  and  Ti  will  discuss  later  in  this  chapter.  They 
are  flow  related  variables. 

The  average  continuity  and  momentum  equations  can  be  linearized.  Noda  et  al.(1974) 
gave  a simple  form  of  those  equations  after  neglecting  the  nonlinear  and  time  dependent 
terms  in  (2.11)  and  (2.20).  They  are 


d 

dxi 


[UiD]  = 0 


(2.23) 


Q=-gD 


drj  IdSij  1_ 

dxi  p dxj  p 


(2.24) 


The  surface  shear  stresses  and  lateral  shear  stresses  are  also  neglected. 


2.2.3  Wave  Energy  Conservation  Equation 

The  equation  considering  both  refraction  and  diffraction  with  current  has  been  devel- 
oped by  Boiij  (1981)  and  corrections  were  made  later  by  Kirby  (1984).  Both  used  varia- 
tional principle.  In  this  study,  a method  that  modified  Kirby’s  derivation  is  developed.  In 
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parallel,  the  perturbation  technique  is  also  applied.  Same  results  are  obtained  and  are  given 
in  Appendix  1.  Here  only  the  derivation  from  variational  principle  is  presented. 

The  variational  principle,  first  applied  in  the  irrotational  flow  motion  by  Luke  (1967), 
states  that  the  variation  of  a certain  conservative  quantity,  often  an  integral,  vanishes  if  the 
function  that  describes  the  evolution  of  the  physical  process  is  subjected  to  small  variations. 
It  is  represented  as 

S[J^J  Ldxdt]  = 0 (2.25) 

where  x it ' + yj  is  in  horizontal  space.  For  irrotational  flow,  the  energy  is  conserved  and 

L can  be  chosen  as  the  Bernoulli  sum: 


(2.26) 


where 

V is  the  three  dimensional  gradient  operator, 

„ d -5  d ^ d . 

t,  j and  k are  unit  vectors  in  x,  y and  z directions,  respectively. 


The  potential  function  <f>  consists  of  two  parts 


<i>  = K + <i>c  (2.27) 

where  <f>yj  and  are  wave  potential  and  current  potential  functions, respectively.  A weak 
point  about  the  expression  of  (j)c  is  that  the  mean  current  field  may  be  rotational,  since 
we  do  not  know  a priori  whether  the  mean  flow  field  derived  from  momentum  equations 
satisfies  the  condition  of  irrotationality  everywhere,  i.e.. 


curl  U = 0 


However,  the  irrotational  assumption  will  be  used  here. 
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For  a slowly  varying  bottom  topography  <f>  can  be  expressed  as  (Kirby,  1983); 

^‘^<f>w2{x,ii''^‘^y,liy,Hx,z,t)  (2.28) 

where  /i  is  a scale  parameter  with  an  order  of  magnitude  0(V/,h/A:/i)  representing  the  spatial 
change  over  the  space  of  wave  length  in  which  V/,  is  horizontal  spatial  gradient  and  A:  = 2;r/L 
is  wave  number  where  L is  wave  length.  8 is  the  strength  parameter  of  current  defined  as 
c is  the  Stokes  wave  steepness  parameter  e = 0{ak),  a is  the  wave  amplitude. 
The  spatial  scale  parameter  fx  is  introduced  to  differentiate  two  different  spatial  scales 
in  the  model,  one  is  related  to  the  fast  wave-like  functions,  and  the  other  one  is  connected 
to  the  slow  spatial  modulations,  such  as  bottom  topography,  the  changes  of  current  speed, 
wave  length  etc.  Therefore,  we  have 

X = x + fix  = X + X2  (2.29) 

y = M^/^y  + My  = -h  ^2  (2.30) 

and  their  derivatives  are  replaced  by 

d d d 

dx  ~ 


dy  ^ dYi 


+ M 


dYi 


(2.32) 


The  scales  as  selected  imply  that  spatial  derivative  with  respect  to  phase  function  occurs 
only  in  the  x-direction  which  is  true  when  the  angle  of  wave  propagating  is  small  with 


respect  to  the  x-axis,  i.e.,  it  is  assumed  that  no  fast  wave-like  variations  occur  in  the  y- 
direction.  It  is  consistent  with  the  forward  scattering  approach  that  the  complex  amplitude 
A is  allowed  to  vary  as  0(/ii/2)  in  the  y-direction.  The  problem  is  further  considered  to  be 
quasi-steady;  time  derivatives  only  appear  in  wave  term. 

There  are  three  magnitude  parameters  8,  /i  and  6 in  the  formulation.  DiflFerent  ratios 
among  those  parameters  will  yield  different  forms  of  the  energy  equation.  In  the  present 
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study,  we  adopt  Kirby’s  (1983)  model  given  the  relationship  as  follows 


S 


1 


(2.33) 


which  represent  slow  varying  topography,  strong  mean  current,  and  small  incident  wave 
angles. 

In  equation  (2.28)  <f>^i  and  <f>yj2  are  first  and  second  order  potential  function,  respec- 
tively. <j>^i  and  <f>^2  are  expressed  as  follows: 


<i>wi  = 

<t>w'i  = A^{x,ti^^^y,nx,ny,t)f^{fix,fj.y,z) 


(2.34) 

(2.35) 


The  can  be  split  into  two  parts:  a modulated  amplitude  function  and  a phase  function: 

A'p{x,tx^l'^y,iix,ny,t)  = A'”"{n^^^y,fj,x,/j.y)e'"^^  + c.c  (2.36) 

where  A'^  is  the  complex  amplitude,  ip  is  the  phase  function  and  c.c  is  the  complex  conju- 
gate. The  phase  function  is  defined  as 

Ip  = j kidx  - (jjt  (2.37) 

where  ki  is  the  wave  number  in  main  (or  x)  direction  and  u = 2nfT  is  the  angular  wave 
frequency  in  which  T is  the  wave  period.  The  real  wave  direction  which  is  different  from 
the  main  direction  is  absorbed  in  the  amplitude  function.  The  subscripts  and  superscripts 
in  the  equations  represent  the  order  and  harmonics,  respectively.  However  the  absolute 
superscript  values  should  be  no  bigger  than  the  subscript  values.  Substituting  (2.34)  and 
(2.35)  into  (2.28),  yields 

<p  = £- Vc(f^x,  e^y)  + eAj"(x,  ey,  e^x,  e^y,  t)f^{e^x,  e^y,  z)  + 

e^A^{x,  ey,  e^x,  e^y,  t)f^{e^x,  e^y,  z)  (2.38) 

where  S and  /x  have  been  replaced  by  and  e^,  respectively,  by  virtue  of  Eq.(2.33). 
Substituting  (2.38)  into  (2.26),  we  obtain  with  (2.37),  term  by  term 


<Pt  = eA^jr+^^A^/r 


(2.39) 
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V<^  = + ef^ATx'i  + e^fr^TyJ  + eVrV^Af  + e^A^V^fr 

+eV2'"A-^;-  + c^fTA^yJ  + e^f^VhA^  + 

+{eA^f^^  + e^A^f^^)k  (2.40) 


where 

V;,  is  the  horizontal  gradient  operator, 


„ d d 
^‘  = ai‘  + a-y’- 


Therefore,  by  definition  we  have 


^h(f>c  = Ui  + Vj 


(2.41) 


where  velocities  include  wave  induced  mean  current. 

The  first  two  terms  in  Eq.(2.40)  represents  linear  wave  superimposed  upon  uniform 
current  over  flat  bottom.  In  Eq.(2.26),  the  (V<^)^  term  when  computed  to  fourth  order 
becomes 


{Vcf>y  = {yh<l>c)^  + €^frf?ATxA^x  + e^frfiATY^A’ly^+e*f^f^ATxA’^jc 

+fVC/i"  + ^^fTzf^zA^A'^  + e2Uf'^A^x  + e^2Vf^A^y^ 

+6^2 f^Vh<l>c  • V/,A;^  + eHA^Vh<t>c  • ^hfr  + + e^2V f^A^y^ 

■ ^hA'^  + e*2A^Vh4>,  • ^hfr  + ^^“^fr/iATxA^x^ 


+^*Vr'f?x,ATA^x  + e^Vrf2ATxA^x  + e^2ATA^f^J^, 


(2.42) 


The  subscripts  X,  Y and  z are  partial  derivatives.  The  superscripts  m and  n represent 
interaction  between  two  harmonics. 

Free  surface  envelope  can  be  considered  as  two  parts:  one  is  related  to  mean  current 
as  r)Q  and  another  one  is  driven  by  waves  as  »/'.  Later  we  will  see  rjo  has  the  same  order  as 

|U|. 


ri{x,y,t)  = rjo{x,y)  + r)'{x,y,t) 


(2.43) 
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or  in  the  stretched  coordinates 

n{x,€y,€^x,e^y,t)  = + n'{x,ey,e^x,e'^y,t)  (2.44) 

where  rf  has  the  order  0(e)  at  most  and  it  can  be  expressed  as 

r?'  = erjr  + e^r?^  + 0(e®)  (2.45) 

where  rji  and  rj2  are  first  and  second  order  surface  elevations,  respectively.  By  virtue  of 
(2.45),  (2.43)  can  be  rewritten  as 


»7  = »?o  + + e^r?^ 


(2.46) 


Expanding  (2.26)  in  Taylor  series  about  z — rjo,  results  in 

i = /j^.  + 

= ly, + + A*, + + 5-?'"  1-1^, + 


+e'''*^l*i  + - ^9A  + --' 


(2.47) 


Introducing  a shifting  coordinate  defined  as 


z'  = z -r)o,  h'  = h + T)o 


and  substituting  into  (2.47),  we  get,  after  dropping  the  primes  on  h'  and  z'  for  convenience. 


- ^gho 


(2.48) 


where  we  denote  ho  as  still  water  depth.  The  expansion  is  carried  out  to  fourth  term  only 
as  the  higher  order  terms  will  not  affect  the  third  order  wave  equation. 
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Certain  integrations  in  (2.48)  that  involve  of  Eqs.(2.34)  and  (2.35)  are  defined  as 
follows  for  future  use: 


r frdz  = Q- 

J — h 

(2.49) 

r frffdz  = 

f 

(2.50) 

' "-h 

(2.51) 

where  t,  j,  m and  n are  free  indexes  and  subscript  z represents  partial  derivatives.  From 
Stokes  wave  theory,  the  well-known  expressions  for  the  fp  and  /I*  are 

fl  = /r^  = 


/I  = f2 


cosh  k{h  + z) 
cosh  kh 

(2.52) 

cosh2A:(h  + z) 
cosh  2kh 

(2.53) 

The  typical  results  in  (2.49)-(2.51)  that  will  be  used  later  are  given  as  follows: 

= - 
9 


pii 

^11 


pl2 

^12 


cosh  kh 
k sinh®  h 

1 + 2 cosh^  kh 
3k  cosh  kh  sinh^  kh 


p22  _ ^ sinh4A:/i, 

-^22  - . , 8 . , (o  + 


sinh®  kh'~2  8k 


R 


11 

11 


- khc,) 


pl2 

-n-12 


4A: 


3 cosh  kh  sinh  kh 
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p22 

^22 


,sinh4A:/i  ,, 
-( h) 


2 sinh®  kn  k 


where  c is  the  wave  celerity,  c = k\ 


Cg  is  the  group  velocity,  Cg  = c • n with 

” ~ + sinh2ifcfc)> 


(T  — y/gktanhkh  is  the  local  angular  wave  frequency,  it  will  be  discussed  in  section  (2.2.4). 
Substituting  (2.39)  and  (2.42)  into  (2.48),with  the  definition  (2.49)-(2.51),  we  have 


L = 


+eQ^UA'^x  + e’^QTyATv,  + e^Q?^h<l>c  • + e^P^^TxAl 


+e  PiT^?xA?Ar  + e^i?iT^r>l2  + + 6^2”"  + + c^r/r/r^u 


+ + e^nTUf'^A'lx  + + ^^Tfi^h<t>c  • V^AJ* 

+eSTUfU2X  + ^STVfiA^y,  + e^j=/r/2"ArxA?;,  + e^.f^J-^ATA^ 

+ \^V2frf^ATxA^,^  + \eSyrj?^ATA^,  + + eSTVf^A’ly^ 

+ ^f"?2”‘/u A?e  + \e\rf^A2t  + f“»7iSr/i"  A?t  + fH^^TxA^x 
+ l^U^fufuzATA-,  + + \e\^V  f^.A'ly  ^ f^,A\x 

A^SXri^VnAlx  + .A?,  + 6'‘if3-£7/[‘,,A?;,  (2.54) 


where  the  condition  of  V/,/”  = 0 has  been  involved.  All  terms  are  evaluated  at  z = 0.  For 
simplicity,  we  only  keep  terms  containing  r}\  and  A^,  since  as  will  be  shown  L varies  with 
respect  to  r}i  and  A\  only.  In  (2.54),  the  superscripts  are  free  indexes  (nonrepeating)  which 
may  not  assume  the  same  value  for  different  terms,  but  they  give  same  value  in  one  term. 
The  terms  of  are  defined  as 

^2  = Smn'nT'li  m + n = k (2.55) 


?3  — SkmnVlVTVl  k + m + n = j 


(2.56) 
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where 

when  m = n or  k = m = n,  Smn  and  Skmn  equal  to  unit,  otherwise  they  assume  the 
values  of  2 and  3,  respectively.  Those  give 


?2  = 


?2  = 


?2^  = 


2riWi 


= 


o„o„i 


?3  ^ ^ + 3rjirjiT]i 


oO^O^-l 


After  introducing  the  following 


+ e^f/2  + e^f/3  + f^f/4  + • • • 


(2.57) 


Equation  (2.54)  can  be  separated  by  grouping  the  terms  of  the  same  order: 


0(e) 


= \qTAu  + QTUaTx  + gmvT  + \rin'^h<t>c? 


(2.58) 


O(e^) 


^2  = \p^^ArxAlx  + + QTVA'^y^  + 

+T)TUf^A’lx 


(2.59) 


O(e^) 


Ls  = QT^h<i>c  • ^hA^  + PTATxA'^.x  + + p»?r»?2”  + riTfU'it 

+ri?fUit  + Iri^ifrf^ATxA’lx  + + r,TVf!^A^y^ 

+vTUf^A^x  + r,TUf^A^,x  + 5?2’”A"  A?,  + h^Uf^.A^^x 


(2.60) 
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There  is  no  Lq  term,  as  it  contains  neither  r}i  nor  term.  For  each  order,  varying 
with  respect  to  rji  and  Ai,  we  will  obtain  a series  of  equations.  The  variational  principle 
gives 

dLi  d ,dLi.  dLi  , 

dB  dt^dBt^ 

where  B represents  either  rji  or  Ai,  Dg  is  the  dissipation  term  that  is  caused  by  bottom  fric- 
tion and  turbulence  that  eventually  are  converted  into  internal  energy,  which  is  irreversible. 
(One  possible  formation  of  is  given  in  Appendix  A by  using  perturbation  method.) 
Equation  (2.62)  can  be  written  in  detail  as 

r_£_  _ d d d d.ddd 

^dB  dX^^dBx,^  = (2.63) 

The  following  discussions  are  based  on  the  orders  of  L. 

1)  0(e) 

Varying  (2.58)  with  respect  to  rjJ”  yields 

»?o  = -^(Vft^c)^  (2.64) 

which  states  that  the  surface  elevation  by  mean  current  is  simply  balanced  by  kinetic  pres- 
sure. The  same  result  is  obtained  if  one  uses  Bernoulli  equation  at  order  0(1). 

2)  0(e2) 

Elquation  (2.59)  varied  with  respect  to  rf'^,  results  in 

»?r  = + 
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= 

9^^  Dt 

^ is  the  total  derivative  related  to  the  x-  direction 


(2.65) 


— = — U— 

Dt~  dX 

For  n = 0 or  the  zeroth  harmonic,  is  independent  of  X and  t,  which  results  in 

77i  = 0 for  n = 0 (2.66) 

From  Stokes  linear  waves,  we  know  that 

(2.67) 

(2-68) 

where  a is  a complex  amplitude  and  ai  is  defined  as 

(Ti=u}-Uki  (2.69) 

which  is  the  apparent  local  angular  frequency  that  represents  the  waves  riding  on  unidirec- 
tional flow.  According  to  linear  wave  theory,  the  dispersion  relationship  is 


a = (j  - Uki-Vk2  (2.70) 

where  k\  kcos9  and  k2  = fcsin^  are  the  wave  numbers  in  x and  y direction,  respectively 
and  6 is  the  incident  wave  angle  as  shown  in  Fig.2.1.  From  (2.69)  and  (2.70),  we  find 


O'!  = (T  + Vk2  = <r-(-0(6) 


since  for  small  incident  wave  angle,  sin0  ~ 0 = 0(^)-  Hence  for  n = 1 and  n = —1,/^^  = 
/i  = 1 at  z = 0,  we  have 


IM  = 

9 Dt  2^ 


(2.71) 
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^ g Dt  2 


(2.72) 


both  are  the  first  harmonic  solutions  and  the  same  as  given  by  linear  wave  theory,  rfi 
and  have  the  same  amplitudes,  but  traveling  in  opposite  directions.  As  mentioned 
before,  the  complex  amplitude  ”a”  absorbs  the  difference  between  real  wave  direction  and 
its  main  direction,  since  phase  function  tp  is  only  defined  in  x the  direction.  Next,  varying 
L2  with  respect  to  A^,  {j  = m or  n),  we  have 

Jffnn  Am  _ pmn  Am  TT  _ n.  in  A,n\ 

"11  /i  »?it  - Ayxx  - ^ h Vix  = 0 (2.73) 

a)  let  m = 0 (the  zeroth  harmonic),  it  gives 

= 0 


Therefore, 

AJ  = 0 or  = 0 from  (2.51) 
Later  we  can  prove  that  if  /i^  = 0,  then  /°  = 0. 
b)  Considering  the  first  harmonic,  we  obtain 

[ccg{kl  - k^)]^  = 0 


(2.74) 


(2.75) 


It  is  important  to  note  that  (2.75)  is  of  the  order  of  O(e^),  since  ^ j^2  ^ j^2  gj^2  g ^ 

0(e  ).  This  term  later  will  be  combined  with  other  third  order  terms  to  form  the  wave 
equation  which  is  of  the  order  of  O(e^) 

3)  O(c^) 

Varying  Z13  with  respect  to  r/J  and  Aj  (j  = n,  m or  k)  respectively,  we  obtain 


+ f2^^t  + \frfiATxA^x  + + V/r  A?y^  + r?r/r.A5*, 


(2.76) 


A”A"  Af  - V* . (QfVk#,)  - F.r  - A“«  - rr/nA^,xnth 
-vf^nTr,  - u/r^Tx  - /uviiS  - A".»f-)5r  = o 


(2.77) 


26 


we  emphasize  again  that  the  superscripts  only  take  the  same  values  in  one  term;  there  is  no 
relationship  among  different  terms.  There  are  only  two  unknowns,  r/j  and  Aj , in  the  above 
two  equations  and  so  that  they  can  be  solved, 
a)  Collecting  all  the  terms  with  zeroth  harmonics  in  (2.76),  we  get 


»?2 


+Vi 


IX -^ix  + fL  ^ + 'll 


ll-^ix 


Here  the  condition  of  f}  = has  been  applied. 

The  values  on  the  right  hand  side  of  above  equation  are  known,  so  that 

0 I , 5A:^tanh^A:h  ktanhkh. 

+ ^ — 1 


(2.78) 


Since  we  only  consider  up  to  O(e^),  fciand  ai  can  be  replaced  by  k and  a,  and  vice  versa. 
Equation  (2.78)  becomes 

This  is  the  set-down  of  mean  free  surface  with  reference  to  the  shifted  coordinates.  From 
(2.77),  we  obtain 

^2  = . (grv^fli,)  (2.80) 

where  m can  assume  any  arbitary  value.  Hence,  A°  can  not  be  determined  as  remains 
unknown.  However  we  will  show  later  that  A§/°^  is  non-zero.  From  (2.80),  the  previous 
statement  in  case  2 is  confirmed,  that  is,  let  m = 0,  we  must  have  /°  = 0 for  = 0. 
Otherwise,  A2  will  approach  infinity. 

b)  The  first  order  harmonics  in  (2.76)  yield  the  following  solution 


'll  — ~ »/2<^1^2] 

From  (2.77),  together  with  (2.81),  we  have 


(2.81) 


^2(^1!' + = 0 


(2.82) 
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For  m — 1,  if  we  eissume  that  /j  has  the  same  form  as  // , the  terms  in  the  bracket  are 

identically  equal  zero.  But,  if  m is  any  value  less  than  or  equal  to  1,  i.e.,  m < | 1 |,  we  may 

select  ^2  = 0 such  that 

n\  = (2-83) 

which  states  that  the  surface  elevation  of  the  first  harmonic  of  the  second  order  of  the  wave 
is  balanced  by  the  variation  of  wave  amplitude  modulation  in  longshore  direction, 
c)  For  the  first  negative  harmonic,  (2.76)  and  (2.77)  give 

= 0 (2.84) 


»72 


-1 

n = — 1 nir  P ^ 


t — ay,  e 
0-1 


(2.85) 


d)  Considering  the  second  harmonic,  after  lower  harmonic  terms  have  been  dropped, 
obtain  from  (2.76)  and  (2.77),  respectively 

9VI  + fiAl  + + f^UAljc  + flr,\A\t  + fLUri\A\x  = 0 


we 


and 


^12 A2  + flz^\A\  - P12A\xx  - 'lit  ~ (Alx'?i)jf  “ ^'llx  ~ flz'li'nit 

-fLnWix  = 0 


Solving  two  equations  for  two  unknowns  A\  and  yield 


■^2 


.3a\  a 


+ o(e^) 
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(2.86) 


2 I a rAjcosh  kh , , « 

(2  + cosh  2A:A)e*^^  + O(e^) 


8 sinh®  kh 


(2.87) 


They  are  the  same  as  the  second  order  Stokes  waves. 
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4)  O(e^) 


Varying  L4  with  respect  to  rj[,  gives 


■ V/.A?  + + frf^ATxA^x  + 

+ri^ifrfi.A^A^x  + n^ifufuzATA^  + + n?U  + \^T  fl^^Al 

+\^TfuzA^iX  + Vf^A'iy^  + r,rV A"  + Vf^A'^y^  = 0 

For  the  first  harmonic,  we  have 

• V;,Ai  + flzflzAi  ^A2  + flzf2zA\A^  + /i/I^LY-^ax  + fu^i^ ' 

n^i  n A~^ 

+/;.»? + A!i,r‘)  + fUL 


{iAlA^'vl  + A\\:')  + + \/Lsl^  = 0 


where 


Dt 


DA~^ 


Dt 


1 3 — 

— = --ae 


The  only  unknown  in  the  equation  is  so  that 


M = -[-^h<f>c  ■ V,(-)  + ■ V^a  + i—\a  |2] 

O’  <Ti  aaci  a 


(2.88) 


in  which 


qk^ 

^ ~ 16sinh^  ~ + 4sinh^  kh  - sinh^  fchcosh^  kh) 


Next,  varying  Z-4  with  respect  to  A\,  also  considering  first  harmonic,  we  get 

flzfh'li^Al  + flJizViA^  + flz{V2A\  + J7|»7f  + flJu{^2A\ 
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+dA ')  + ‘>l2  + + dni ')  - PnMy,y, 

-(P,Utx)x,  - Pii‘/iUa-  - V*  ■ - fi(,,^‘Alj,)x  - (n°2A\x 

~'l\^ix)x  — flz(^2-^ix  + ?2-^ix)  “ ^ ^2^1  — 0 

where  again  fl  = = 1 is  applied.  Each  term  in  the  above  equation  has  been  calculated 

before.  After  some  tedious  algebraic  manipulation,  we  obtain 

2{cc,ki  + aiU)ax,  + 2cr,VaY,  - i{cc,  - V^)ay,Y,  + 

V 

+ (— )k2]o  + ikccgk'l  a\'^a  — 0 (2.89) 

where  k'  — k^  — + 8-2  tanh^  kh 

Cg  Ssinh'^kh 

ki  and  <ti  in  (2.89)  can  be  replaced  by  k and  a without  introducing  errors  larger  than  O(e^). 
Equation  (2.89)  contains  information  of  wave  refraction  and  diffraction,  as  well  as  wave 
amplitude  dispersion  and  wave-current  interaction.  Finally,  as  we  have  noted  before  the 
first  harmonic  for  Li  varying  with  respect  to  Ai  as  expressed  by  (2.75)  is  third  order  and, 
therefore,  should  be  combined  with  (2.89)  to  reach  the  wave  euqation.  Hence,  we  obtain  in 
the  (x,y)  coordinates 

+ ^,U)a.  + 2^^Va,  - i{cc,  - V^)a„  + + (H) 

<Ti  <Ti 

+tccg{kl  - A:^)]a  + ikccgk'\  a |^a  = 0 (2.90) 

This  is  the  final  form  of  the  wave  equation  without  dissipation.  If  dissipation  is  to  be 
considered,  then  a dissipation  term  De  could  be  included  in  the  left  hand  side.  One  possible 
formation  of  De  is  given  in  Appendix  1.  Following  Kirby  and  Dairy mple  (1983),  a phase- 
shifted  wave  amplitude  is  introduced  to  simplify  the  wave  number  computation 

where  ko  is  the  reference  wave  number.  Equation  (2.91)  implies  that 

A"  = 


(2.92) 
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where  a"  is  the  real  wave  amplitude.  The  final  wave  equation,  after  dropping  the  double 
primes,  takes  the  following  form: 

2{ccgh  + (TiU)A^  + 2aiVAy  - i(cc,  - V^)Ayy  + 

+ {y)y  + iccg[k\  - e)  + 2i{ccgk^  + axU){ko  - A:i)]A  + i{al  - a^)A 
(i  — l)gk‘^  I,  u , , 

where  the  cubic  amplitude  term  in  Eq.(2.90)  has  been  replaced  by  t(a^  - a^)A  (Kirby  and 
Dalrymple,(1986)),  in  which 

Oq  = gk{l  + fie^Dk)  tanh{kh  + /2e) 

it  will  be  discussed  in  next  section.  The  last  term  on  the  left  hand  side  is  the  decay  of 
energy  due  to  wave  breaking  and  the  details  of  which  will  be  discussed  in  section  2.2.6,  and 
W is  defined  as 

{W  = 0 before  wave  breaking 

2^2 

^ “ 4;^)  3.fter  wave  breaking,  7 = 0.39  and  q — 0.4 

The  dissipation  term  De  is  given  in  appendix  A.  Equation  (2.93)  is  the  final  form  of  the 
wave  equation.  It  is  similar  to  the  wave  equation  given  by  Kirby  (1984)  except  the  ad- 
ditional imaginary  term  in  the  second  last  term,  (see  Liu  1986).  For  plane  beach  and  in 
the  absence  of  current,  linearized  (2.90)  reduces  to  the  conventional  one-dimensional  wave 
energy  conservation  equation.  It  reads 

a = O’o{cgo/cg){ho/h) 

where  quantities  with  subscripts  0 refer  to  known  condition  and  h is  the  width  between 
adjacent  wave  orthogonals. 

2.2.4  Dispersion  Relation 

When  waves  travel  over  a varying  topography  or,  are  embraced  by  current,  the  wave 
number  vector  will  change  its  value.  The  four  equations  established  in  the  previous  section 
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enable  us  to  solve  the  four  unknowns:  rf,  U,  V and  H , provided  the  wave  number  vector  is 
known.  Since  wave  number  is  a horizontal  vector,  in  general  two  equations  are  required  to 
solve  its  magnitude  and  direction.  The  wave  number  irrotationality  and  the  wave  disper- 
sion relationship  derived  from  dynamic  free  surface  condition  are  commonly  used.  In  the 
presence  of  a current  field  the  mathematics  could  become  complicated  if  the  nonlinear  be- 
havior is  fully  accounted  for.  In  the  present  study,  a semi-nonlinear  dispersion  relationship 
is  considered  and  is  incorporated  directly  in  the  wave  equation;  thus  the  linear  dispersion 
relationship  is  used  in  the  model.  Also,  by  assuming  small  incident  wave  angle,  the  di- 
rectional information  on  wave  number  is  incorporated  into  the  complex  wave  amplitude  as 
presented  in  the  previous  section,  thus,  avoiding  the  inclusion  of  the  wave  number  irrota- 
tional  equation  explicitly.  The  linear  dispersion  relation  with  current  is  given  as  (Dean  and 
Dalrymple,  1984) 

u = (T  + Uki+Vk2  (2-94) 

where 

(T  - v^^A:tanh  kh 
ki  = k cos  6 
k2  = k sin  6 

In  shallow  water,  tanh  kh  — ► kh  approximately,  Eq.(2.94)  becomes 

u)  = ky/^+Uki  + Vk2  (2.95) 

Several  modifications  have  been  proposed.  Whitham(1967)  gave  a Stokes  nonlinear  disper- 
sion relation  for  intermediate  depth  region  in  the  absence  of  current: 

(To  = ayJl  + e^Dk  (2.96) 

where 


€ = ak 

cosh  4kh  + 8-2  tanh^  kh 

Uk  = 


8 sinh^  kh 
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As  waves  enter  shallow  water,  (2.96)  becomes  invalid  ais  approaches  (kh)~* 

In  the  shallow  water  region.  Hedges  (1976)  proposed  a different  formula 

(^'o  = '^gkta.nhk{h  + a)  (2.97) 

for  the  case  of  no  current  in  which  a is  wave  amplitude.  For  the  same  wave  frequency,  it 
yields  longer  wave  length  than  that  given  by  the  linear  dispersion  relation. 

More  recently,  by  combining  both  cases  discussed  above,  Kirby  and  Dalrymple  (1986) 
constructed  a composite  model  which  satisfies  both  intermediate  water  depth,  Eq.(2.96), 
and  shallow  water  depth,  Eq.(2.97), 

(^"o  = [gk{^  + he^Dk)  tanh(A:/i  + /2f)]^^^  (2.98) 

where 

fi  = tanh®  kh 
fi  = (A:h/sinh  khY 

In  deep  water  region,  (2.98)  approaches  (2.96)  as  — » 0;  conversely,  in  shallow  water 

region,  it  becomes  (2.97)  a.s  /i  — » 0.  Now  assuming  there  is  no  nonlinear  current  effect  on 
wave  dispersion,  the  presence  of  current  causes  a Doppler  shift  and  Eq.(2.98)  becomes 

= ff"o  + Uki  + Vk2  (2.99) 

Figures  2.2  and  2.3  shows  the  results  of  w for  a constant  non-dimensional  velocity  U /cq  = 
—0.1  and  varying  of  e from  0.01  to  0.4,  where  cq  = >/g/k  is  the  deep  water  wave  celerity. 
The  dispersion  equations  given  by  (2.94),  (2.96)  and  (2.97)  are  also  plotted  in  Figs.2.2  and 
2.3  for  comparison. 

The  emprical  nonlinear  dispersion  (2.98)  can  be  used  to  modify  the  wave  equation  as  we 
did  in  Eq.(2.93),  in  which  the  Eq.  (2.98)  is  incorporated  directly  in  the  governing  equation 
in  the  form  of  a term  which  directly  produces  the  required  amplitude  dispersion  (Kirby  and 
Dalrymple,  1986).  Consequently,  the  linear  dispersion  relation  (2.94)  will  be  used  in  this 


model. 
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Figure  2.2:  Variations  of  Dispersion  Relation  with  U jc^  = —0.1  and  e,  - 
Hedges;  - - - Witham; Composite 
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figure  2.3;  Variations  of  Dispersion  Relation  with  UjcQ  — —0.1  and 
Hedges; Whitham; Composite 


Linear;  — 
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Once  the  dispersion  equation  (2.99)  is  brought  into  the  model,  the  cubic  amplitude  term 
in  Eq.(2.90)  should  be  dropped,  since  the  amplitude  dispersion  has  already  been  includeed 
in  the  nonlinear  dispersion  equation  (2.98). 

A wave  number  vector  which  is  irrotational  can  be  defined  as: 

k = Vw  (2.100) 

where  k is  a wave  number  vector,  k = {A:i,A:2}  and  u corresponds  to  the  scaler  phase 
function.  Using  the  mathematical  property  that  the  curl  of  a gradient  is  identically  zero,  it 
is  shown  that 

V X Vw  = 0 (2.101) 

or 

V X k = 0 (2.102) 

In  general,  Eq.(2.102)  together  with  Eq.(2.99)  is  necessary  to  solve  the  wave  number 
field  of  a domain. 

However,  in  the  present  study  the  wave  angle  information  has  already  been  included  in 
the  complex  wave  amplitude  given  by  Eq.(2.92).  Therefore  equation  (2.102)  does  not  need 
to  be  included  explicitly.  However,  the  assumption  that  wave  approaching  angle  is  small 
must  be  observed. 

2-2.5  Radiation  Stresses,  Lateral  Mixing  Stresses  and  Bottom  Stresses 

In  the  wave  momentum  equations  (2.21)  and  (2.22),  there  are  some  undetermined  terms: 
radiation  stresses,  lateral  mixing  stresses  and  bottom  shear  stresses.  They  are  not  indepen- 
dent variables  and  are  calculated  by  the  following  equations. 

Radiation  Stresses 

The  radiation  stresses  are  given  in  (2.19): 

= f pSijdz  + f uf  ujdz  — -pgD^Sij  - piiiUjD 

J — h J — h ^ 
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(2.103) 


The  form  was  simplified  by  Longuet-Higgins  and  Stewart  (1962)  to  second  order  for  a linear 
progressive  wave  system.  For  normal  incident  waves,  it  can  be  approximately  by 

E{2n  - 0.5)  0 

0 E{n  - 0.5) 

where  E is  local  wave  energy  E = ^pg\  a |^,  o is  complex  wave  amplitude  and  n is  the  ratio 
of  wave  group  velocity  to  the  wave  celerity. 

In  general,  when  waves  enter  beach  with  an  angle,  5,y  can  be  obtained  by 


■ s;. 

s.“,  ■ 

:5.y]  = [A][5?.][A] 


where  [ A ] is  directional  matrix. 


A = 


cos  & — sin  6 

sin  9 cos  9 


hence 


Sij]  = 


Sxx 

Sxy 

F?[n(cos^  0 + 1)  - |] 

^En  sin  20 

Sxy 

■ 

\En  sin  20 

£J[n(sin^  0 + 1)  - |] 

where  9 is  local  wave  angle. 

Bottom  Shear  Stresses 

The  bottom  shear  stresses  are  customarily  expressed  as 


(2.104) 


(2.105) 


(2.106) 


m = pfuti  I Uf  I (2.107) 

where  Ut  is  total  velocity  composing  of  both  wave  orbital  and  mean  current  velocities,  and 
uti  is  its  component  form.  The  quantity  / is  the  friction  factor.  The  magnitude  of  total 
velocity  | ut  | is  equal  to  \Zu^  + The  u and  v velocity  components  can  be  expressed  as 

u = U + u'^  = U + u^  cos  9 (2.108) 


v^V  + ul  = U + u'^sin9 


(2.109) 
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where  U and  V are  mean  current  components  as  defined  before,  u™  and  are  wave  orbital 
velocities  in  the  x and  y direction,  respectively  and  = [(u“)^  + 

The  total  velocity  | ut  | can  be  expressed  as 

I Ut  1=  [U^  + V^  + {u'^f  + 2Uu'^  cos9  + 2Vu'^  (2.110) 

The  local  wave  orbital  velocity  u’"  is  purely  oscillating  in  time 


u = Um  cos  at 


(2.111) 


where  Um  is  the  maximum  wave  orbital  velocity  at  the  bottom 


u™,  = 


a a 


sinh  kh 

The  time  averaged  bottom  shear  stresses  in  the  x and  y directions  then  become 

_ J 


/ fT 

’’’bx  = J {u  + UmCOsO  cosat)  I Ut 


dt 


(2.112) 


(2.113) 


— f r 

% = io  ^ I 

Alternative,  for  computational  convenience,  it  can  be  rewritten  as 

f /•2ir 


dt 


Hx  = P-^  {U  + UmCosO  cosat)  I Ut  1 d(at) 


(2.114) 


(2.115) 


f 

^by  = P—  {U  + Umsinecosat)  \ ut  | d{at)  (2.116) 

These  integrals  can  be  calculated  by  applying  Simpson’s  rule.  The  serious  problem  by 
applied  the  above  bottom  friction  model  is  the  computational  time.  A simplified  model  is 
introduced  with 


m = pf  I u't  I Ui 


(2.117) 


I u't  \=  ViP+V^  + Urr, 


where 
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since  there  is  no  interaction  term  between  wave  orbital  velocity  and  mean  current  velocity. 
In  the  component  form  the  above  equation  gives 

= pf\u't\U  (2.118) 


ny  = pf  \ u't\V 


(2.119) 


The  difference  between  the  two  models  (2.107)  and  (2.117)  is  minor  for  cases  U ■>  Um  as 
shown  by  the  example  in  Fig. 2. 4,  in  which  the  same  bottom  friction  coefficient  is  used  (the 
testing  conditions  refer  to  Fig.5.3  in  chapter  5),  but  the  saving  in  computational  time  is 
more  than  60%.  In  fact  the  small  difference  can  be  offset  by  increasing  the  bottom  friction 
coefficient  in  (2.117)  to  account  for  the  absence  of  term.  Equation  (2.117)  will  be  used 
in  our  model. 

Another  very  important  factor  is  selection  of  the  friction  coefficient  /.  Based  on  field  ob- 
servations and  laboratory  experiments  Longuet- Higgins  (1970)  suggested  a value  about  0.01 
for  long-shore  current  calculations.  The  bottom  friction  coefficient  for  waves  alone  presents 
a more  difficult  situation.  For  smooth  turbulent  flow  over  fixed  bed,  Kamphuis(1975)  gave 


/ = 


1 

We 


(2.120) 


where  Rg  is  Reynolds  number,  Rg  — and  u is  the  kinematic  viscosity. 

For  rough  turbulent  flow,  friction  coefficient  no  longer  depends  upon  Reynolds  number. 
Investigations  were  conducted  by  Jonsson(1966)and  Grant  & Madsen(1979).  The  most 
commonly  used  relation  is 


f iexp[5.213(^)0  i9'‘  - 5.977]  ^ < 0.63 

f = ) (2.121) 

[ 0.15  5 > 0.63 

where  r is  the  roughness  of  the  bed  and  a'  =|  a |.  Equation  (2.121)  is  an  approximate 
explicit  solution  given  by  Swart  (1974)  to  an  implicit,  semi-empiricaJ  formula  by  Jonsson 
and  Carlsen  (1976). 
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Figure  2.4:  Velocity  Profiles  for  different  friction  models  — - (2.107);  - - - (2.117) 
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The  situation  becomes  further  complicated  once  combined  wave  with  current,  the  state 
of  the  art  of  modeling  such  flows  is  rather  rudimentary.  The  total  velocity  in  (2.107)  becomes 
much  larger  than  the  wave  orbital  velocity  in  the  present  case  of  very  strong  current.  The 
actual  values  of  / used  in  the  model  will  be  adjusted  as  a function  of  (a,  T)  based  on 
experimental  results  to  be  discussed  in  chapter  4. 

In  application,  the  bottom  friction  coefficient  may  be  treated  differently  in  x and  y 
directions,  in  the  present  study,  the  bottom  friction  coefficient  in  y direction  is  taken  as  one 
half  of  the  value  in  x (or  main)  direction. 

Lateral  Mixing 

Lateral  mixing  due  to  turbulent  diffusion  is  a complicated  flow  phenomenon.  Rigorous 
mathematical  treament  is  beyond  the  state  of  art.  Figure  2.5  shows  the  effects  of  lat- 
eral mixing  on  longshore  current  distribution  (Longuet-Higgins,  1970  a.b).  The  common 
assumption  is  that  turbulent  mixing  is  proportional  to  the  local  mean  velocity  gradient.  Fol- 
lowing Longuet-Higgins’  formulation  on  longshore  current  mixing,  the  lateral  shear  stress 
is  written  as 

dU  dV 

where  and  are  the  horizontal  eddy  viscosities  in  the  x and  y direction  respectively.  We 
define 

= ^wi  + i = x,y  (2.123) 

where  e^i  and  are  the  eddy  viscosities  of  wave  and  current,  respectively  , which  are 
assumed  to  be  independent.  The  inclusion  of  Cd  is  based  on  the  fact  that  even  though 
there  is  no  wave,  the  flow  from  the  inlet  still  generates  turbulence.  As  the  simplest  possible 
approach,  Longuet-Higgins  proposed  a model  for  that  is  proportional  to  the  distance 
offshore,  x,  multiplied  by  velocity  scale 


^WX  = Nxy/^ 


(2.124) 
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Figure  2.5:  Longuet-Higgen’s  Longshore  Current  Profiles 
where  N ia  & dimensionless  constant  whose  limits  Longuet-Higgins  gave  as 


0 < iV  < 0.016 


(2.125) 


Following  Ebersole  euid  DiJrymple  (1979),  N is  chosen  as  0.01  and  is  allowed  to  vary 
hnearly  with  x to  the  breaking  line.  From  this  point  seaward  the  coefficient  remeiins  at 
this  constant  value.  There  is  no  elaborate  consideration  of  c«,y.  Ebersole  and  Dairy mple 
assumed  a constant  value  amd  here  it  is  chosen  as  half  the  value  of  e^x  throughout  the  field. 

The  eddy  viscosity  theoretically,  will  change  from  point  to  point.  One  of  the  formulas 
weis  due  to  Von  Karman’s  similarity  hypothesis: 


(2.126) 


zcy 


= »^2  I 


dy 


(2.127) 
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where 


and  A:  = 0.4  is  Von  Karman’s  constant.  In  the  numerical  calculation,  these  coefficients  go 
to  infinit  near  a point  of  inflection.  An  upper-limit  value  should  be  imposed.  Alternatively, 


both  eddy  viscosities  can  be  treated  as  constants.  Figure  2.6  give  the  effects  of  eddy  viscosi- 
ties on  the  velocity  profiles  (the  testing  conditions  refer  to  Fig. 5. 3 in  chapter  5).  It  shows 
that  they  have  only  minor  difference  between  constant  0.001  rri^/sec  and  Von  Kormon’s 
mixing  model  with  upper-limit  equal  to  0.001  rn^jsec.  Therefore  in  this  investigation,  in- 
stead of  Eqs.(2.126)  and  (2.127)  the  constant  eddy  viscosity  is  selected. 

2.2.6  Wave  Breaking  Criterion 

When  waves  enter  shallow  water,  they  will  eventually  become  unstable  and  break, 
dissipating  the  energy  in  the  form  of  turbulence  and  work  against  bottom  friction.  For  a 
very  mildly  sloping  beach,  spilling  breaker  occurs  and  the  simplest  criterion  predicted  by 
solitary  wave  theory  is  (McGowan  (1894)) 


where  D is  total  water  depth  and  subscript  h denotes  the  value  at  breaking. 

Dally  et  al.  (1984)  proposed  that  the  decay  of  energy  flux  with  distance  in  the  surf 
zone  is  proportional  to  the  excess  of  energy  flux  above  a stable  value,  given  the  relation 


where  q is  a.  constant  related  to  the  rate  of  energy  decay.  The  quantity  [Ecg)g  is  the  stable 
energy  flux  for  a broken  wave  in  water  of  total  depth  D.  The  last  term  in  the  wave  equation 


(2.128) 


(2.129) 


(2.93)  which  appears  only  after  wave  breaking 


2{ccgki  + aiU)W  A 
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Figure  2.6:  Variations  of  Eddy  Vicosity  Coefficients.  — - = 0.  ; •••  =0.001;  =0.01; 

Von  Karman’s  Model  with  upper-limit  =0.001.  unit=  m*/s 


can  now  be  completely  defined  as 
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Q , 

by  virtue  of  Eqs.  (2.128)  and  (2.129).  The  value  of  q is  chosen  as  0.4  following  Kirby  and 
Dalrymple  (1986). 


2.3  Summary 

In  this  chapter,  the  analytical  model  of  a current-wave-topographic  interaction  system 
is  formulated  under  the  assumption  of  strong  current  and  slowly  varying  bottom.  Basic 
equations  together  with  essential  terms  are  developed.  They  are  summarized  as  follows: 
Depth-Integrated  Continuity  Equation 


^ dUD  dVD  _ 
dt  ^ dx  ^ dy  ^ 


(2.131) 


Depth-Integrated  Momentum  Equation  in  ^-Direction 


— + U—  +V—  = -n^  - _ J_- 

dt  dx  dy  ^dx  pD  dx  pD  dy  ^ p dy  pD^’’^ 

Depth-Integrated  Momentum  Equation  in  y-Direction 


^ + U—  + V—^-n^-±^-±^  + h^_±- 

dt  dx  dy  ^ dy  pD  dx  pD  dy  ^ p dx  pD^^'^ 
in  which  we  have  substituted  (2.131)  into  (2.21)  and  (2.22),  where 


cO 

^xy 

£^[n(cos^  0 + \)  - \] 

\En  sin  20 

oO 

^yy  . 

\En  sin  26 

E[n(sin^  0 -h  1)  - |] 

dU 


dV 


(2.132) 


(2.133) 


(2.134) 

(2.135) 


— Ctut  + €(.,•;  is  a constant  and  e^x  is  given  by 
^wx  = Nx\/gD  with  N = 0.01; 
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Wave  Equation 


2(cc,fci  + aiU)A,  + 2aiVAy  - i(cc,  - V^)Ayy  + 

V 

+ (~)v  + — k^)  + 2i{ccgki  + aiU){ko  — fci)]A  + »(ctq  — ct^)A 

(»  — l)pA:^  / 1/ 

+ + = “ (2138) 

where 


= [5^(1  + /ie^£>*.)tanh(A:/i  + /2e)]^^^ 

/i  = tanh®  kh 

f2  = {kh/  sinh  kh)* 

cosh  4kh  + 8-2  tanh  kh 

” : — j““ 

8 sinh^  kh 

= 0 before  wave  breaJcing 

^ ="  ^(1  “ after  wave  breaking,  7 = 0.39  and  q = 0.4  (2139) 

When  waves  enter  the  surf  zone,  the  last  term  in  (2.138)  will  be  brought  into  calculation 
till  waves  reach  stable  condition. 

The  dispersion  relation  used  in  this  model  is 

uj  = <T  + Uki +Vk2  (2.140) 

and  the  wave  angle  is  included  in  the  complex  wave  amplitude  A. 


CHAPTER  3 
NUMERICAL  SCHEME 


In  this  chapter,  the  numerical  schemes  for  solving  the  continuity,  momentum  and  wave 
equations  are  presented,  along  with  boundary  and  initial  conditions. 

The  numerical  scheme  requires  (1)  the  definition  of  a grid  system  with  a system  which 
provides  a systematic  method  for  identifying  variables  of  interest  and  (2)  the  conversion  of 
the  governing  equations  into  finite  difference  forms. 

3.1  Grid  System  and  Definition  of  Variables 

A rectangular  grid  mesh  was  established  over  the  area  of  interest  as  shown  in  Fig. 3.1, 
where  x and  y denote  the  onshore  and  longshore  direction,  respectively.  Two  kinds  of 
variables  are  employed  in  the  computation:  the  grid-center  variable  and  the  grid-side  vari- 
able. Practically  all  the  quantities  of  interest  are  represented  by  the  grid-center  form  which 
expresses  the  value  at  the  center  of  the  grid  and  has  the  following  notation: 

where  X is  the  grid-center  variable,  the  subscripts  i,  j represent  the  ith  grid  side  in  the 
x-direction  and  the  jth  grid  side  in  the  y-direction,  respectively,  and  the  superscript  n 
represents  the  nth  time  step.  The  grid-side  variable  expresses  the  value  of  the  quantity 
along  any  grid  side.  The  only  variable  that  appeared  in  this  form  in  the  present  study  is  the 
velocity  vector;  it  is  denoted  by  a subscript  s to  differentiate  it  from  the  grid-center  value 
and  has  the  following  form: 
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Ys  is  the  vector  quantity  along  a grid  side:  the  subscripts  are  such  that  t denotes  the  »th 
grid  point  in  the  x-direction  on  the  left  of  the  vector  and  j denotes  the  jth  grid  point  in 
the  y-direction  on  the  right  of  the  vector. 

In  the  present  study,  the  grid-side  velocity  vector  was  applied  in  the  continuity  and 
momentum  equations  and  a grid-center  velocity  vector  was  applied  in  the  wave  equation. 
The  definitions  of  grid-side  and  grid-center  variables  are  shown  in  Fig. 3. 2 

3.2  Finite  Difference  Scheme  for  Combined  Continuity-  Momentum  Equation 

The  Crank-Nicolson  finite  difference  scheme  in  double  sweep  form  is  applied  here  to 
solve  the  continuity  and  momentum  equations  and  later  also  to  solve  the  wave  equation. 

3.2.1  Linearized  Implicit  Model 

A variety  of  numerical  schemes  are  presently  available  to  solve  the  equations  for  nearshore 
hydrodynamics  similar  to  those  developed  in  Chapter  2.  Ebersole  and  Dalrymple  (1979) 
presented  an  explicit  model  for  nearshore  circulation  which  requires  small  time  steps  due 
to  the  stability  criterion,  thus  relatively  long  computational  time.  Vemulakonda  (1982) 
presented  an  implicit  model  allowing  for  much  larger  time  steps.  Both  models  adopt  a 
three  time  level  schemes.  Butler  (1981)  and  Sheng(l981)proposed  a two  time  level  implicit 
model.  This  method  not  only  saves  computational  time  but  also  shortens  the  development 
of  the  difference  scheme  for  the  governing  equations.  To  demonstrate  the  structure  of  their 
model,  the  simplified  linearized  equations  are  discussed  here.  The  depth  averaged  linearized 
continuity  and  momentum  equations  are 


dr)  ^ dQ^  , dQy 

T 1 — — == 


dt  dx  dy 


(3.1) 


(3.2) 
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Figure  3.1:  Grid  Scheme 
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where  Q^.  = UD  and  Qy  = VD. 
matrix  form,  linearized  equations 


dt 


U and  V in  this  section  denote 
can  be  rewritten  as 


(3.3) 


the  grid-side  values.  In 


Wt  + AW^  + BWy  = 0 


(3.4) 


where 


( ^ 

( ° 

1 

0 ^ 

( ® 

0 

1 A 

w = 

Qx 

A=\  gD 
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0 B = 

0 

0 

0 

< Qy 

J 

1 0 

0 

oj 

0 

o; 

Sheng  (1981)  implemented  a two  time  level  fully  implicit  scheme  in  computing  the  external 
model  in  his  three  dimensional  mode. 


'M  It,] 


AT^  + - {Qy)m  = 0 (3.5) 


At 


(3.6) 


iQy)i,r  - 

At 


+ 


_ ^+1  1 


= 0 


or  in  the  matrix  form 


(3.7) 


+ + = 0 (3.8) 

where  6x  and  6y  are  central  spatial  difference  operators  which  provide  the  differentiated 
quantities  at  the  same  points  as  the  time  difference  quantities.  Equation(3.8)  involves  only 
the  (n-l-1)  time  level  in  the  last  two  terms,  i.e.,  the  whole  equation  is  staggered  in  time. 
Rearranging  (3.8)  by  time  levels,  we  obtain 


(1  + 2/3x  + 2/3y)W’^+^  = ly" 


(3.9) 


„ A — -;n, — B6^ 


2Ax 


2Ay 


where 
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By  adding  the  quantity  of  — W'^),  which  has  a second  order  truncation  error 

as 


we  have 


{MfAB 


d^W 

dtdxdy 


(1  + 2^,)(1  + = (1  + 4My)W^  (3.10) 

This  factorized  finite  difference  equation  is  still  a second  order  approximation  to  the  differ- 
ential equation  (3.4).  The  advantage  of  using  (3.10)  lies  in  the  fact  that  the  solution  of  the 
factorized  form  can  be  split  into  two  separate  one  dimensional  operations.  To  achieve  this, 
we  introduce  an  intermediate  value,  W*,  such  that 


{1  + 2^^)W*  = {l-20y)W 


(3.11) 


(1  + 2/3„)W""+^  = W*  + 2/3yW 


(3.12) 


From  now  on,  quantities  without  the  time  level  denoted  are  the  values  at  time  level  n. 
Equation  (3.10)  is  recoverable  from  (3.11)  and  (3.12)  by  eliminating  W*. 

A double  sweep  method  is  used  to  solve  (3.11)  and  (3.12).  At  the  first  sweep,  the  values 
at  time  level  n are  known  we  solve  for  IF*.  The  full  solution  then  is  obtained  based 

upon  known  W*  values.  After  some  manipulations,  the  first  sweep  in  x direction  gives 


~ Qx  + ^^gOSxTj*  = 0 


At 


and  the  second  sweep  in  the  y direction  yields 


(3.13) 


(3.14) 


I ^ - ^SyQy  = 0 

-Qy  + ^QDSyr^-^  = 0 

On  the  first  sweep,  and  rf*  are  found,  and  on  next  sweep  the  values  of  and 

are  obtained.  Therefore,  after  one  loop,  the  quantities  at  time  level  n -M  are  known.  The 
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advantages  of  this  scheme  are  (a)  the  finite  difference  equations  are  brief,  at  second  sweep 
(3.14)  there  is  no  Qx  involved;  (b)  a significant  savings  in  computation  time  is  realized  with 
large  time  steps  and  less  computations. 


3.2.2  Numerical  Scheme  for  the  Full  Equations 


Following  the  above  model,  the  full  equations  (2.131),  (2.132)  and  (2.133)  in  chapter  2 
can  be  written  in  numerical  form  as  (see  Appendix  C) 
x-sweep 


(1  + AtCj;)f7”'*'^  — U + Atrrix  = 0 

y-sweep 


(3.15) 


^+1  -rf*  + - ^fi{VD)  = 0 


(3.16) 


(1  + Atc„)F"+i  - y + ^gSyfj^+^  + Atrrii  = 0 
where  c^U  and  CyV  are  bottom  shear  stresses  in  x and  y direction,  respectively.  By  (2.136) 
in  chapter  2,  we  have 


= and  cy  = ^\u't\y 

where  | U(  |,-  is  the  total  velocity  defined  in  (2.137)  and  here  is  specified  in  each  direction, 
mi  and  m2  are  the  remaining  terms  in  Eqs.  (2.132)  and  (2.133). 

Again,  on  the  x sweep  the  solutions  for  ff*  and  [7"+^  are  obtained;  these  up-dated 
values  are  brought  into  the  y sweep  to  obtain  and  V^+\  Writing  the  above  two  set 
of  equations  in  detail,  we  have 
a)  X-Sweep 
Continuity  equation 


- ^ij)  + (A+i.;  + A,)  - + A_i.y)] 

2 Ay  ~ + A,j-i)]  = 0 


(3.17) 
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X-Direction  Momentum  Equation 

-tmr'  - v,j) + 


+ l V,j  ,)  li-lj) 

p{Dij  - A-l.j)  ~ i^xy)i,j-l  + (Sj^),_i  y+1  - (S'iy),_iy_i] 

~ ( £}.  J -|_  Di-I  j)  ^ + i^^yx)i,j  (3.18) 


8Ay 
1 


where 


I (“<)».J  \x  ~ [(“wi)t,j  + (’^m)t-l,;]  + 2(Ui)l,  J 

Um  is  wave  orbital  velocity  as  defined  in  (1.112)  in  chapter  2. 

Kki  = {t^/!y+[0.25(V.-y+V.-y+i  + V;_i.y  + V._i,y+0]'}'^^ 


{^^yy)i,i  - 4(Ay)2'f(^‘-i+i  ~ ^»jl((f«'v)<,y  + (^««v),-,y+i  + ktt,i/),_i  y + kiuv),_i,y+i) 

+2((Ccj/),- y + (fey),-  y^i)]  ~ {Ui  j - l^t,y-l)[((etuy),_y  + ktuy),_y_i 

+ (ft«v),_l,y  + (ei«y),_i_y_i)  + 2((ecy),-,y  + (fey),- y_i)]} 


(FLyx) 


^-i,i+i)[((^t«i)f,j+i  + {^wx)ij  + (c«;z),_iy 
+ 2((e„)..y^i  + (e..),._i^.+i)]  - (V.- y - V.--i,y)[((c,,).  . 

+ (fiui),_lj  + (eiuz),_y_i  + (f  wx  )t-ij-i)  + 2((eci),_y  + (eci),_iy)]} 


b)  Y-Sweep 
Continuity  Equation 
1 


~ 2 Ay  ~ ^■,y(A’,y  + A-,j_i)*]  = 0 (3.19) 

Y-Direction  Momentum  Equation 

A(  '''•■>  + 5aJ ^ 
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(Vt.j+1  Vij  i)  Vi, til) 

-(•^a^J.-lj  + (5x„),+  l,j-l  - (5jj,),_l,y_l]  - --  ^ ^ ^ ^^Syy)iJ 


p{^i,j  ”1"  •^»’,j— i)*  ^y 
iSvv)i.j-i)  - 2(D  . + D - + {FI'xx)i,j  (3.20) 

V ' *»J  1/ * 


where  the  bottom  shear  stress  coffficient  has  taken  as  half  the  value  of  that  in  x the  direction 
as  stated  in  (2.2.5),  and 


I (“<)«, J ly  — [(*^m)i,;  + («m)i,j-l]  + 2{uy)l,j 

MiJ  = {Kj  + [0.25(C//71  + C^/;+i.  + ^^^^^^ 

{FLxx)t,j  — {(^+l.J  ~ ^■,i[((fiui),+ij  + (Ciuz),- y + (etuz),+  i_y_i  + (^’"*)«,j-l) 

+ 2((ecz),-|-i,y  + (fcz)j,y)]  ~ (Vi.y  - Vj-i  j)[((e^j;),- y + (eujj:),._j  y 
+ (etuz)j,y_i  + (fwz),_i,y_i)  + 2((eci),.y  + (fez), _i,y)]} 


{FLxy)i,j 


+ Mi,j-i)  + 2((ec„),+i,y  + (^o„),+i,y_i)]  - {u::/'  - 
+ (^’"t/).,>-l  + (fiuv),_l,y  + (fiuy),_i_y_i)  + 2((ecy),.y  + (ecy),j_i)]} 


In  (3.17)  to  (3.20),  the  quantities  with  no  time  level  are  at  time  level  n;  and  the  subscript 
(♦)  outside  the  bracket  in(3.19)  and  (3.20)  denotes  up-dated  values,  i.e.,  Z)*  = h + ff*.  Also 
we  have  used  the  latest  x direction  velocity,  C/"+i,  in  the  y sweep  to  accelerate  the  conver- 
gence. In  application,  the  {ccx)  and  (f^y)  are  considered  as  constants  as  dicussed  in  chapter  2. 


3.2.3  Method  of  Calculation  and  Boundary  Conditions 

A double  sweep  method  has  been  developed  in  each  direction,i.e.,  in  both  x-sweep  and 
y-sweep.  The  advantage  of  this  sweep  method  lies  in  the  fact  that  it  gives  an  explicit  formula 
at  each  sweep  to  calculate  the  results. 
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Collecting  the  terms  from  (3.17)  to  (3.20)  by  time  level,  so  that  the  unknown  quantities 
are  to  the  left  and  the  known  qantities  are  to  the  right,  we  have 
X- Sweep 


Y-Sweep 


- PI  v!j  = B<.i 


(3.21) 


(3.22) 


where  the  quatities  without  superscripts  are  at  time  level  n.  The  definitions  of  undefined 
quantities  are  given  in  Appendix  D,  which  also  shows  the  details  of  the  following  derivation. 
The  X sweep  equations,  after  some  manipulations,  can  be  reorganized  as 

= Fli  + 


Ur^J-^  = F2i_i  + 


(3.23) 


where 


Fh 

Eli 

F2i 

E2i 


Fj,jF2i 
1 + Pi,jE2i 
Fj-i,j 
1 + PijE2i 
Ej+i,}  ~ 

A+ij  + 

R1 

A+ij  + RlEli^i 


The  boundary  conditions  state  that  at  the  shoreline  (i=idry),  the  velocities  should  be  equal 
to  given  velocities  in  the  inlet  section  and  they  disappear  elsewhere,  and  the  mean  water 
displacement  r?,  assumes  to  be  identically  equal  zero.  Therefore,  the  first  sweep  will  solve 
coefficients  FI,  F2,  El  and  E2  starting  from  i=idry  till  offshore  line  (i=l).  Once  those 
coefficients  are  determined,  the  second  sweep  will  directly  yield  results  of  rf*  and  17"+^ 
starting  from  i=l  (offshore  boundary)  till  i=iwet,  where  iwet  = idry-1,  is  the  last  grid  at 
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which  water  depth  has  non-zero  positive  value.  The  offshore  boundary  condition  is  given  at 
deep  water,  i=l,  by  specifying  the  mean  water  displacement  to  be  zero.  This  is  approximate 
assumption  to  maintain  a steady  state  mean  flow  in  the  region  of  concern  with  no  increase 
or  decrease  of  water  mass.  A more  rigorous  condition  is  to  allow  for  water  level  variation 
at  the  offshore  boundary  such  that  the  change  in  hydrostatic  force  is  balanced  by  the  net 
momentum  flux  and  bottom  shears  in  the  x-direction. 


The  lateral  boundary  conditions  basically  have  three  kinds  for  the  concerned  region: 
the  symmetric  free  boundary  conditions;  the  no-flow  boundary  conditions  and  the  non- 
symmetric  free  boundary  conditions.  They  will  be  specified  for  specific  problems  later. 

The  y sweep  is  similar  to  the  x sweep  and  (3.22)  can  be  rewritten  as 


where 


(3.24) 


'*.7 

= F3j  + 

F4j-i  + 

FZj 

_ Eij  — QijF4j 

1 + QijE4j 

EZj 

_ 

1 + Qi,jE4j 

F4j 

_ “ R2F3j+i 

Rij+1  -f  R2E3j^i 

E4j 

R2 

Ri,j+i  "b  R2E3j^i 

The  first  sweep  will  determine  F3,  E3,  F4  and  E4,  then  the  following  sweep  will  give  the 
values  of  and  The  shoreline  boundary  condition  for  V is  set  equal  to  zero  at 

i=idry;  and  offshore  boundary  condition  at  i=l  supposes  that  it  satisfies  the  mass  con- 
servation, i.e.,  {VD)x  = 0.  The  lateral  boundary  conditions  will  be  discussed  for  different 


situations  later. 
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3.2.4  Stability  Criteria 

The  stability  criteria  of  these  schemes  can  not  be  established  theoretically.  In  general, 
the  stability  criterion  for  explicit  scheme  can  be  guided  by  the  following  inequality: 


< V(Ax)"  + (Ay)^ 

I I +VqD 

where  | ut  | is  the  magnitude  of  total  velocity.  For  application,  it  becomes 

At  < 

^ V9D  ^ 


(3.25) 


(3.26) 


for  I Ut 

Actually,  the  critical  time  step.  At,  in  our  study  by  using  a two  time  level  implicit  model 
is  at  least  one  order  larger  than  that  given  by  (3.26). 


3-3  Finite  Difference  Scheme  for  the  Wave  Equation 

The  final  form  of  wave  equation  (2.138)  in  chapter  2 is  mixed  parabolic-hyperbolic 
partial  dififerential  equation.  The  advantage  of  parabolic  type  partial  differential  equation 
IS  that  the  marching  method  can  be  used  in  step  by  step  calculation.  Booij  (1981)  proposed 
a parabolic  model  following  on  Radder  (1979).  Then  the  wave  field  can  be  split  into  a 
transmitted  and  a reflected  field.  In  his  parabolic  approach,  the  reflected  field  is  completely 
neglected.  Kirby  (1983)  gave  another  approximation,  after  assuming  that  the  complex 
amplitude  can  absorb  the  difference  between  real  and  image  wave  directions  for  small  angle 
incident  waves,  and  also  assuming  that  wave  amplitudes  are  changing  more  rapidly  in  y 
direction  than  in  x direction  for  slowly  varying  topography.  Both  of  those  assumptions 
have  been  introduced  in  our  derivation  of  the  wave  equation. 

For  the  parabolic  partial  differential  equation,  the  Crank-Nicolson  finite  difference 
scheme  is  applied.  For  simplicity,  (2.138)  in  chapter  2 can  be  rewritten  as 

A.  + CoA,  - iCiA,,  + + (C4  + iC5)A  = 0 (3.27) 
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where 


in  which 


Co 

Cx 

C2 

C3 

C4 

Cs 


2axV 


ccg  - 


V_ 

0^1 

R + W 


— ~ + (^0  - ki)-  R + -^(o-Q  - <T^) 


P = 2{ccgki  + aiU) 

R = [Y^ 

P cosh^  kh  V 2(jj 


As  the  waves  enter  the  surf  zone,  the  second  term  in  C4  will  be  added  into  calculation.  The 
quantity  of  | A in  is  the  value  evaluated  at  last  iteration,  unlike  Kirby  who  specified 
this  value  at  the  same  iteration.  The  values  of  velocity  components  U and  V are  grid-center 
values  in  this  calculation. 


The  Crank-Nicolson  finite  difference  scheme  constructs  Eq.(3.27)  as 

~ Ai  ’ ^ - A+lj-i)  + (Co)i,y 

(^■j+i  Aj-i)]  - 2(Aj/)2  [(^i)»+i,i('^+ij+i  ~ ^^+i,j  + A,+i  j_i)  -1-  (Cl),- j 


(A...  - 2^, + + a<,) 

A,-+i,y  + 


2Ax  {{C2)i+ij  + (C2),-j) 

^((^3)«-fl,j+l  - (C3),-+i  y_i)  ^ ^ ((C3),-J+1  - (Cs),- j_i) 


(C'2).-.y 


4Ay^  (C^2).+i,i 

+ 2 [((^4).+l,i  + t(C5),-+ij)A,-+i_y  + {{C4)ij  + t'(C5),-_y)A,-j-]  = 0 


'4«,yJ 


(3.28) 


A similar  double  sweep  method  is  used  to  solve  (3.28),  resulting  in 


Aij  — i?2j_i  + iJlj-jA,- y_l 


(3.29) 


The  definitions  of  Rl,  R2  and  other  details  are  also  given  in  Appendix  C. 
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The  boundary  conditions  are  stated  as  follows: 

a)  at  offshore,  i=l,  the  wave  amplitudes  are  specified  in  complex  forms; 

b)  at  shoreline,  i=idry,  the  wave  amplitudes  are  supposed  to  die  off  completely; 

c)  the  lateral  boundary  conditions  can  be  considered  as 


Ay  = ik  sin  OA 


(3.30) 


which  assumes  that  the  waves  travel  mainly  in  the  x direction  on  a slowly  varying  topography 
and  Snell’s  law  is  valid  for  first  order  correction.  The  finite  difference  form  of  (3.30)  gives 


— Ai,j 


1 + Q 

1-g 


where 


(3.31) 


g = (j'A:sin^^),j+i 


The  local  wave  propagation  angle  9,  as  discussed  in  chapter  2,  can  be  found  through 
the  complex  wave  amplitude.  By  Eq.(2.92)  of  chapter  2,  we  have 

(3-32) 


= «2A:2Ay  (3.33) 

where  the  double  primes  in  Eq.(2.92)  have  been  dropped.  The  local  wave  angles  are  in- 
cluded both  in  A:i  and  k^,  so  that  either  of  above  equations  can  be  used  to  obtain  6. 


3.4  Method  of  Solution 

With  the  specified  initial  and  boundary  conditions,  the  finite  difference  equations  (3.21), 
(3.22)  and  (3.28)  can  be  solved  in  a loop  of  iteration. 

The  initial  condition  is  assumed  to  be  the  state  of  rest.  The  velocity  field,  both  U and  V 
components  as  well  as  the  mean  free  surface  displacement,  fj,  are  initially  set  equal  to  zero.  A 
bottom  topography  is  input  into  the  model.  The  initial  wave  amplitude  field  is  calculated 
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by  solving  the  wave  energy  equation  with  velocities  equal  to  zero,  which  corresponds  to 
Snell’s  law  for  plane  beach  before  waves  break. 

To  avoid  "shock  loading”  the  model,  both  wave  amplitude  and  velocities  of  inlet  dis- 
charge are  built  up  from  zero  to  their  termination  value  over  a number  of  iterations.  We 
choose  100  iterations  to  set  up  those  values  by  using  the  hyperbolic  tangent  function 

F = Fq  tanh(/T/25)  (3.34) 

where  F is  either  wave  amplitude  or  velocity  of  inlet  discharge;  Fq  is  their  termination 
values,  IT  is  the  number  of  iteration.  The  approach  declares  that  the  termination  value 
will  be  reached  as  total  iteration  is  up  to  100,  while  tanh4  approaches  the  value  of  unity. 

The  program  will  run  till  the  variables  U,  V,  rj  and  A reach  their  steady  state.  The 
convergency  is  also  checked  by  the  balance  of  toted  discharge  in  and  out  of  the  controlled 
boundary. 

There  are  two  coefficients  to  be  determined  in  a running  case:  one  is  the  bottom 
friction  coefficient  and  the  other  is  the  eddy  viscosity  coefficient.  For  fixed  inlet  discharge, 
the  eddy  viscosity  is  assumed  to  be  constant,  but  the  bottom  friction  coefficient  may  vary 
for  different  combinations  of  wave  and  current.  Further  explanations  on  how  to  select  these 
two  coefficients  are  given  in  the  next  two  chapters. 


CHAPTER  4 

PERFORMANCE  OF  THE  MODEL 


To  evaluate  the  performance  of  the  numerical  model,  a number  of  tests  were  run,  mainly 
utilizing  existing  numerical  models  of  simple  cases  as  bench  marks.  Little  attempt  was  made 
to  compare  directly  with  existing  experimental  data  for  two  reasons: 

(1)  Experimental  data  that  can  be  used  to  test  the  model  to  its  full  term,  i.e.,  current- 
wave  interaction  due  to  inlet  on  sloping  beach,  does  not  exist. 

(2)  The  bench  mark  numerical  models  selected  here  usually  already  did  their  compar- 
isons with  laboratory  data,  whenever  available.  Since  the  purpose  of  the  present  model  is 
not  to  seek  improvement  of  any  specific  numerical  model  but  rather  to  extend  the  capa- 
bilities to  new  applications,  it  is  deemed  suflFcient  for  verification  purpose  to  compare  with 
numerical  results  only. 

4.1  Wave  Setup  on  Plane  Beach 

The  model  was  run  for  a case  of  normal  incident  wave  on  a plane  smooth  beach.  The 
conditions  were  as  follows:  T=1.14  sec,  deep  water  wave  height  Hq  = 6.45  cm  and  beach 
slope  = 0.082.  This  was  the  same  case  used  by  Ebserole  and  Dalrymple  (1979)  to  compare 
with  laboratory  data  by  Bowen,  Inman  and  Simmons  (1968).  Since  it  is  a two  dimensional 
case,  we  have  selected  m = 52  in  the  x-  direction  (on-offshore)  with  Aa:  = 0.25m  and  n = 4 
in  the  y-direction  (longshore)  with  Ay  = 0.4m.  The  water  depth  at  the  deep  water  end 
in  the  numerical  model  was  about  1 m and  the  mean  water  level  displacement  due  to  set 
down  can  be  neglected.  The  results  are  shown  in  Fig.4.1.  The  wave  height  was  computed 
on  the  basis  of  linear  wave  shoaling  and  the  wave  equation  as  expressed  by  (2.138)  was 
used  for  wave  height  computation.  The  data  from  Bowen,  Inman  and  Simmons  (1968)  were 
also  plotted  for  comparisons.  Apparently,  these  two  approaches  are  similar  except  around 
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Figure  4.1:  Set-up  on  a Plane  Beach.  Eq.(2.138); Linear  Wave;  A Data 
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breaking  line. 

^ Open  Lateral  Boundary  Conditions  and  Longshore  Current  Distributions  on  a Plane  Beach 

Two  typical  open  boundary  conditions  can  be  specified  for  the  present  numerical  model. 

One  of  them  is  the  periodic  lateral  boundary  conditions,  which  requires  that 

Q(«-,i)  = g(f,AT  + i) 
g(i,2)  = g(,-,AT  + 2) 

Q{i,S)  = Q{i,N  + 3) 

and  so  forth,  where  Q is  any  quantity  of  interest.  This  boundary  condition  can  obviously 
be  applied  to  topographies  with  periodicity  such  as  the  example  given  by  Ebsersole  and 
Dalrymple  (1979),  but  can  also  be  relaxed  to  apply  to  an  arbitrary  topography  that  has  no 
periodicity  as  long  as  the  boundaries  are  chosen  far  enough  away  from  the  area  of  interest 
and  there  is  no  source  or  sink  in  the  area. 

In  the  case  where  a source  or  sink  exists,  such  as  the  present  study  with  an  inlet, 
the  periodic  lateral  boundary  condition  can  no  longer  be  applied  as  the  out-flow  from 
downstream  would  not  be  balanced  by  the  in-flow  from  the  upstream.  An  approximate  free 
boundary  condition  in  this  case  is  to  let 

dy 

at  far  enough  distance  on  both  upstream  and  downstream.  The  above  boundary  condition 
implies  no  longshore  variations  far  away  from  the  area  of  interest. 

A typical  example  for  open-boundary-condition  computation  is  the  longshore  current 
generated  by  obliquely  incident  waves.  The  problem  was  studied  analytically  by  Longuet- 
Higgins  (1970).  The  case  illustrated  here  is  for  the  conditions:  beach  slope  = 0.075;  incident 
wave  angle  = 5°;  bottom  friction  coefficient  = 0.01;  wave  height  = 0.02  m and  wave  period 
= 1 second.  The  results  are  shown  in  Fig. 4. 2. 
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If  there  is  no  lateral  mixing  effects  at  all,  the  velocity  decreases  rapidly  from  maximum 
to  zero  at  breakline.  The  apparent  existence  of  velocities  outside  the  breaker  zone  is  due 
to  the  fact  that  the  advective  acceleration  terms  in  the  differenced  y momentum  equation 
cause  velocities  outside  the  surf  zone  to  be  contaminated  by  velocities  inside  the  breaking 
zone.  The  strong  mixing  effects  are  demonstrated  in  the  other  two  cases:  one  only  has  the 
lateral  mixing  due  to  breaking  with  a mixing  coefficient  = 0.01  m^/sec  and  the  other  in- 
cludes both  mixing  from  breaking  waves  with  coefficient=  0.01  m^/sec  and  flow  turbulence 
with  an  eddy  viscosity  coefficient  = 0.005  m*/sec. 

4.3  Nearshore  Circulation 

In  this  example,  currents  induced  by  waves  over  irregular  bottom  topography  are  com- 
puted. The  bottom  topography  is  slightly  modified  from  that  used  by  Noda  et  al.  (1974) 
and  Ebersole  and  Dalrymple  (1979)  and  is  generated  by  the  following  equation: 

h = s ■ a:[  1 -H  sin^°  ^(y  - >/itan^)]  (4.1) 

where 

s:  beach  slope  = 0.025 
x:  distance  from  shoreline 
A:  periodic  beach  length  = 104  m 
A:  amplitude  of  bottom  variation  = 10  m 
= 30° 

in  which  ^/x  is  used  instead  of  x on  the  right  hand  side  as  in  the  original  equation  by  Noda. 
This  modification  is  necessary  to  reduce  the  variations  on  the  lateral  boundary  conditions. 
We  choose  m=40  in  x-direction  with  Ax  = 5.0  m and  n=27  in  y-direction  with  Ay  = 4.0 
m.  The  other  parameters  are  given  as 

T = 4.0sec.;  oq  = 0.75m;  ^ = 5°;  / = 0.012 

breaking  mixing  coefficient  = O.Olm^/sec. 
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Figure  4.2:  Lonshore  Current  Generated  by  Obliquely  Incident  Waves.  — no  mixing; 

breaking  mixing  coeff.=;0.01mV«;  ' - " breaking  mixing  coefF  =0.01m*/s  and  eddy  vis. 
coefF. =0 .005 /« 
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The  resulting  velocity  field  is  shown  in  Fig. 4. 3.  Although  this  case  offers  no  direct  com- 
parison with  the  result  by  Ebersole  and  Dalrymple,  the  flow  patterns  are  similar  and  the 
current  field  is  physically  plausible. 


This  example  deals  with  the  wave  refraction  problem  as  illustrated  in  the  previous 
three  cases,  but  also  the  combined  effects  of  refraction  and  diffraction  under  the  influence 
of  topography  and  the  interaction  with  current.  A case  of  a submerged  circular  shoal  with 
a parabolic  configuration  resting  on  a beach  is  studied.  Similar  topography  is  commonly 
seen  around  inlet  where  an  ebb  shoal  is  located  outside  an  inlet.  This  idealized  topography 
was  first  investigated  by  Radder  (1979)  and  then  by  Kirby  and  Dalrymple  (1983)  on  a flat 
bottom  with  normal  incident  waves.  Therefore,  the  numerical  model  was  first  run  at  the 
same  condition  used  by  Radder.  Then,  the  case  of  obliquely  incident  waves  was  considered 
with  the  shoal  on  a plane  beach.  Finally,  the  model  was  treated  with  the  presence  of  an 
inlet  to  the  above  configuration. 

4.4.1  The  Circular  Shoal  on  a Flat  Bottom 

Radder  (1979)  suggested  two  configurations  of  shoals  in  his  study  of  parabolic  wave 
equation.  Configuration  I,  of  a circular  shoal  on  a flat  bottom  is  represented  by  the  depth 


4.4  Combined  Refraction,  Diffraction  and  Current  Interaction  Case 


profile 


r < R 
r>  R 


(4.2) 


where 


r^  = {x-  Xm)^  + {y- 
eo  = {ho  - hm)/R^ 


Xm  and  ym  are  the  coordinates  of  the  center  of  the  shoal,  and  is  the  water  depth  at  that 
point;  ho  is  the  water  depth  at  the  flat  bottom  area.  The  portion  of  r < R prescribes  a 
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Figure  4.3:  Velocity  Field  of  Periodic  Bottom  Topography 
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parabolic  curve.  The  numerical  values  selected  are 


hm/i?  = 0.0625;  ho/i?  = 0.1875;  Lo/R  = Q.h 


where  Lq  is  the  wave  length  in  the  flat  bottom  area. 

For  such  a configuration,  the  wave  rays  of  a linear  refraction  approximation  are  focus 
into  cusped  caustics.  Peregrine  (1983)  and  Kirby  and  Dalrymple  (1983)  discussed  this 
case  and  have  shown  that  two  jumps  in  wave  condition  fan  out  in  approximately  the  same 
manner  as  the  caustics  emanating  from  the  region  of  the  cusp.  The  region  between  the 
jump  conditions,  directly  in  the  lee  of  the  submerged  shoal,  would  then  be  dominated  by 
waves  of  approximately  uniform  amplitude  propagating  in  the  incident  wave  direction. 

The  grid  spacings  used  by  Radder  are  described  as 


and  it  is  noted  here  that  the  larger  the  value  of  Lq/Ax,  the  better  the  resolution.  For  the 
shoal  studied  here,  Lq  = 0.25m  is  used,  the  corresponding  Z/o/Ay  = 2.5.  The  center  of  the 
shoal,  with  the  radius  R = 10Ax(5Ay)  is  located  along  the  line  of  symmetry  with  respect 
to  the  y-axis. 

Without  consideration  of  current,  the  wave  energy  equation  becomes 

2CCgklAx  - iCCgAyy  + [ {klCCg)^  + icCj(A:J  - k^)  + i2{kiccg){ko  - ki)  ]A 


Following  Kirby  and  Dalrymple  (1983),  the  non-dimensional  incident  wave  amplitudes  are 
chosen  as  fcoAo  = 0.16  and  0.32  in  which  ko  = 2k / Lq.  The  variations  of  wave  eimplitudes 
along  the  center  line  {y/R  = 0)  in  the  x direction  are  shown  in  Fig.4.4  in  which  x/R  = 0 
starts  from  the  fore-edge  of  the  shoal  with  the  center  of  shoal  located  at  x/i?  = 1.0.  For  both 
cases,  in  the  neighborhood  of  the  center  of  the  shoal  the  wave  amplitudes  drop  significantly 
below  the  incident  wave  amplitude.  The  peak  value  of  koAa  - 0.32  is  less  than  half  the  value 
of  koAo  = 0.16;  it  indicates  the  effects  of  nonlinearity  with  the  increase  of  wave  steepness. 


Ay/Ax  = 0.5;  and  Z-o/Ax  = 1, 2, 4 and  8. 


(4.3) 


R/flO 
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Figure  4.4:  Wave  Amplitude  Variations  along  Centerline Aoito  = 0.16; Aoko  = 0.32 
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In  other  words,  the  wave  energy  is  much  easier  to  be  trapped  and  harder  to  decay  for  the 
shorter  waves  at  the  line  of  symmetry.  Also,  there  is  a slight  phase  shift  in  the  position  of 
maximum  wave  amplitudes.  Figures  4.5  and  4.6  gives  the  wave  amplitude  variations  in  the 
longshore  direction  at  the  location  of  x/i?  = 3.0,  4.0,  5.0  and  6.0,  respectively.  y/R  = 0is 
the  symmetric  center  line  of  the  shoal.  These  plots  show  that  as  waves  travel  downstream 
from  the  shoal  the  influence  zone  becomes  wider.  The  spreading  of  the  influence  zone 
for  both  cases  are  almost  the  same.  Since  wave  amplitudes  for  ifcoAo  = 0.16  have  a large 
ampliflcation  a.t  y/R  = 0,  which  can  also  be  shown  in  Fig.4.4,  they  drop  rapidly  to  below 
the  incident  wave  amplitude  to  keep  the  wave  energy  balance. 

These  results  are  simileu-  to  those  described  by  Kirby  and  Dalrymple  (1983). 

4.4.2  Circular  Symmetric  Shoal  on  Plane  Beach 

In  the  case  of  the  circular  shoal  located  at  the  center  of  a uniform  sloping  beach  as  shown 
m Fig.4.7,  the  water  depth,  referring  to  Fig.4.8,  is  prescribed  by  the  following  formulas; 

{h  = X tan  a r > Rcos  a 

h = xtana-h'  r < Rcosa 

where 


r^  = (x-X„)2  + (y-y„)2, 


h'  = 


h"  = 
B = 
C = 


h" / cos  a, 

-B  + VB^  - 4C 


2 

R? 

2r 

hn tan^  a 

sin  a 

r2 

R^ 

sin^  a tan^  a ’ 


and 


The  values  used  in  the  examples  are:  R = 0.7521m;  = 0.085m;  Ax  = Ay  = 0.15m  and 

the  center  of  shoal  is  located  at  grid  point  i=17  and  j=13. 

The  model  is  first  run  for  a case  with  no  inlet  and  then  for  a case  with  the  presence  of 
an  inlet.  In  both  cases,  the  input  waves  form  an  angle  of  5°  with  the  normal  of  the  shoreline. 


OH/b  " ' Ob/b 
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Figure  4.5:  Wave  Amplitude  Variations  across  Centerline 


Ao^o  = 0.16; AoA:o  = 0.32 
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Figure  4.6:  Wave  Amplitude  Variations  across  Centerline Aoko  = 0.16; AoAtq  = 0.32 
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Figure  4.7;  Bottom  Topography  : a Circular  Shcal  with  Parabolic  Configuration  on  the 
Plane  Beach 
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Figure  4.8:  Description  of  Water  Depth 
The  other  input  conditions  are 

r = 1.0sec.;  oo  = 0.02m;  a = 0.075;  / = 0.01 

breaking  mixing  coefficient  = 0.01m* /sec. 

For  the  case  of  no  inlet,  the  eddy  viscosity  coefficient  is  equal  to  zero  and  for  the  case  with 
inlet,  a value  of  0.005m*/sec  is  used.  The  inlet  discharge  is  equal  to  0.007m^/scc. 

The  velocity  field  with  no  inlet  is  shown  in  Fig.4.9.  Waves  break  at  the  top  of  the 
shoal.  The  current  induced  by  the  waves  is  characterized  by  a strong  onshore  flow  over  the 
shoal  consistent  with  the  wave  incident  angle;  at  the  lee  side  of  the  shoal,  the  flow  splits 
into  two  directions,  one  upstream  and  one  downstream.  The  upstream  flow  develops  into 
two  vortexes  rotating  in  opposite  directions;  the  downstream  flow  forms  a longshore  flow 

system  with  two  portions.  This  longshore  current  pattern  at  the  downstream  boundary  is 
sketched  in  Fig. 4. 10. 

The  velocity  field  with  the  inlet  is  shown  in  Fig. 4. 11.  The  effects  of  inlet  flo%v  are  seen 
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to  be  quite  pronounced.  First  of  all,  the  two  upstream  vortices  are  strengthened  and  pushed 
farther  offshore.  Secondly,  the  downstream  longshore  current  system  is  also  strengthened 
with  a flow  reversal  close  to  the  shoreline,  which  is  caused  by  the  lower  mean  water  surface 
elevation  at  the  mid-section  of  the  shoreline  in  the  shadow  of  the  shoal  and  also  caused  by 
the  lower  pressure  due  to  the  higher  velocities  in  the  inlet,  and  the  influence  is  extended 
further  offshore  which  is  due  to  the  fact  that  the  momentum  of  the  inlet  discharge  pushes 
the  flow  into  the  lower  pressure  region.  The  mean  water  surface  elevation  is  also  lowered  in 
this  region.  The  longshore  current  distribution  at  the  downstream  boundary  for  the  case 
with  inlet  is  also  shown  in  Fig.4.10.  The  location  of  null  flow  on  the  lee  side  of  the  shoal  is 
now  being  pushed  seaward  onto  the  shoal  due  to  the  inlet  flow  momentum. 

Figures  4.12  and  4.13  show  the  contour  lines  of  wave  amplitudes.  The  unit  in  these 
figures  is  10  * meter.  For  the  case  of  no  inlet,  the  wave  refraction  effect  causes  the  waves 
to  focus  on  the  shoal  and  this  leads  to  wave  breaking  and  a low  region  just  leeward  of  the 
shoal  crest.  The  diffraction  efiect,  on  the  other  hand  causes  waves  to  focus  on  the  leeward 
of  the  shoal  and  creates  a local  high  region  there. 

The  effects  of  current  are  most  pronounced  in  the  center  region  between  the  inlet  and  the 
shoal.  The  local  high  region  behind  the  shoal  as  mentioned  above  is  enhanced  considerably. 
The  waves  at  the  location  of  inlet  entrance  are  also  steepened  and  the  heights  increased. 
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Figure  4.9:  Velocity  Vector  Field  for  Circular  Shoal  Located 


on  the  Plane  Be2ich  (no  inlet) 
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Figure  4.10:  Longshore  Current  Distributions  at  Lowest  Boundary.  No  Inlet;  - - - with 

Inlet 
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Figure  4.11:  Velocity  Vector  field  for  Circular  Shoal  Located  on  the  Plane  Beach  with 
Current  from  Inlet 
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Figure  4.12:  Contour  Plot  of  Wave  Amplitudes  on  the  Circular  Shoal  Located  on  the  Plane 
Beach  (no  inlet) 
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Figure  4.13:  Contour  Plot  of  Wave  Amplitudes  on  the  Circular  Shoal  Located  on  the  Plane 
Beach  with  Current  from  Inlet 


CHAPTER  5 

EXPERIMENTAL  COMPARISONS 


In  order  to  test  the  validity  of  the  numerical  model,  laboratory  experiments  were  con- 
ducted in  a square  basin  in  the  laboratory  of  the  Coastal  and  Oceanographic  Department, 
University  of  Florida.  The  experimental  set  up,  procedures  and  results  together  with  nu- 
merical simulations  are  reported  in  this  chapter. 

5.1  Purpose  of  Experiments 

Although  there  are  m2iny  theoretical  studied  in  recent  years  on  the  subject  of  wave- 
current  interaction,  laboratory  experiments  and  field  data  axe  scarce.  Ismail- Awadalla 
(1981)  conducted  laboratory  experiments  in  wave-current  interaction  on  a flat  bottom  with 
a nozzle  injecting  dye  against  the  normal  incident  waves.  The  results  demonstrated  the 
physical  phenomenon  of  current-wave  interaction  but  were  not  sufficiently  quantitative  to 
verify  numerical  models. 

In  the  present  experiments,  we  focus  on  a beach-inlet  system,  where  the  phenomena  of 
wave  transformation  over  a gentle  slope,  coupled  with  the  effects  of  a non-uniform  current 
from  an  inlet  can  be  studied.  There  are  four  basic  quantities  to  be  investigated:  mean 
velocities:  Z7  in  x direction  and  U in  y direction;  the  mean  water  surface  displacements 
T)  and  the  wave  amplitudes  A.  Since  ^ is  a small  quantity  and  requires  a long  time  to 
establish  a steady  value,  the  effects  of  wave  reflection  in  the  wave  basin  preclude  meaningful 
measurements.  Only  wave  amplitude  and  current  field  were  investigated. 
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5.2  Facilities  and  Instruments 
5.2.1  The  Laboratory  Facilities 

The  experiments  were  conducted  in  a square  7.22  m by  7.22  m wave  basin  with  a depth 
of  0.8  m as  shown  in  Fig. 5.1  and  5.2.  A fixed  plane  beach  was  constructed  of  cement  with 
1:13.3  slope  (or  s = 0.075).  The  toe  of  the  beach  was  0.85  m from  the  wave  generator.  The 
inlet  was  located  at  the  center  of  the  beach  end.  It  was  1.44  m wide  and  the  bottom  was 
horizontal.  The  total  length  from  the  leading  edge  of  the  inlet  to  the  diffuser  was  1.2  m 
as  shown  in  Fig.5.2.  For  all  cases  tested  the  water  depth  was  kept  0.369  m measured  from 
the  bottom  of  the  basin.  Thus,  the  water  depth  was  0.045  m in  the  inlet.  A grid  system 
was  established  as  shown  in  Fig. 5.3  to  facilitate  numerical  computation  and  laboratory 
measurement.  The  grid  size  is  0.12  m x 0.18  m.  There  are  43  grids  in  the  x-direction  and 
40  in  y-direction.  The  inlet  was  located  at  the  center  of  basin  from  j = 17  to  j = 24.  The 
leading  edge  of  the  inlet  was  at  i = 37.  The  toe  of  the  sloping  beach  was  located  at  i = 1. 

Waves  were  generated  by  a flap  type  wave  generator.  The  bottom  of  the  flap  was 
elevated  0.2  m above  the  bottom  to  permit  the  excess  water  from  the  inlet  to  return  freely 
under  it.  The  wave  generator  was  driven  by  a D.C.  motor.  Wave  period  and  amplitude 
were  controlled  by  adjusting  the  frequency  of  the  motor  and  the  eccentricity  of  the  linkage. 
Uniformly  long  crested  waves  can  be  obtained  for  small  amplitude  waves.  Cross  waves 
would  appear  as  the  wave  amplitudes  became  larger.  Standing  waves  eiIso  existed  under 
certain  circumstance  due  to  the  short  distance  from  the  toe  of  the  beach  to  the  wave  flap. 

To  prevent  the  leakage  of  wave  energy  from  behind  the  wavemaker,  horse  hair  mats  were 
placed  between  the  wave  paddle  and  the  wall,  so  that  the  wave  motion  behind  the  paddle 
was  greatly  reduced.  On  the  other  hand,  due  to  space  limitation,  wave  energy  absorber 
could  not  be  placed  at  the  beach  end.  Therefore,  reflection  of  waves  was  serious  at  times. 

The  current  system  consisted  of  a pipe  system,  a flow  filter,  a diffuser  and  a reservoir. 
As  shown  in  Fig.5.1,  a 2.5  horsepower  centrifugal  pump  located  between  the  flow  filter  and 
the  diffuser  was  used  to  transport  the  water.  The  flow  filter  situated  behind  the  wave  flap 
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Figure  5.1;  Plaji  View  of  Wave  Tank 
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Figure  5.2:  Cross-Section  Views  of  Wave  Tank 
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Figure  5.3:  Cross-Sections  of  Meeisurement 
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has  ten  control  valves,  allowing  the  water  to  flow  naturally  to  the  recirculation  pipe  without 
disturbing  natural  flow  condition.  A diffuser  with  a length  1.4  m was  installed  in  a trench 
at  the  end  of  the  inlet.  It  was  a section  of  pipe  with  a series  of  different  sized  holes  to 
distribute  the  flow  uniformly.  The  gap  around  the  diffuser  was  filled  in  with  pebbles  and  a 
filter  screen  made  by  two  layer  a net  filter  with  horse  hair  in  between  was  set  up  in  front 
of  the  trench  to  promote  uniform  flow  across  the  inlet.  However,  the  distance  between  the 
difipuser  and  the  intersection  of  the  shoreline  and  the  inlet  entrance  was  only  0.6  m.  It  was 
very  difficult  to  attain  a uniform  flow  condition  for  such  a short  distance.  After  many  fine 
tunings,  an  acceptable  uniform  flow  condition  was  finally  reached.  Fig. 5.4  shows  the  velocity 
distribution  in  the  inlet  at  the  cross-section  that  intersects  the  shoreline.  The  adopted  final 
velocity  distribution  was  7/25/86  profile. 

5.2.2  Instrumentation 

The  main  properties  measured  included  wave  characteristics  and  current.  Four  resis- 
tance type  wave  gages  were  used  to  measured  the  spatial  distributions  of  wave  heights  and 
periods.  Current  was  measured  by  a Marsh-Mcburney  type  electro-magnetic  current  meter. 
In  addition  to  wave  and  current,  the  discharge  from  the  inlet  was  monitored  by  a commer- 
cial flow  meter  manufactured  by  Sitnet  Scientific.  The  wave  gage  was  calibrated  statically 
by  varying  the  gage  submergence  stepwise.  It  was  found  that  the  gage  could  not  maintain 
linearity  for  the  full  gage  length.  Therefore,  gage  calibration  was  finally  done  for  each  lo- 
cation prior  to  each  test.  The  current  meter  was  supplied  with  calibration  but  was  further 
calibrated,  verified  and  adjusted  by  towing  test  in  a long  flume  in  the  same  laboratory. 
The  flow  meter  was  also  calibrated  by  measuring  the  discharge  filling  in  a reservoir.  Fig.5.5 
shows  the  calibration  curves  of  the  three  instruments;  the  results  were  consistent  and  stable. 

A 6-channel  strip  chart  recorder  was  used  to  display  the  data.  The  wave  gage  signal  was 
first  amplified  before  feeding  into  recorder.  The  current  meter  signal  was  time-integrated 
to  eliminate  short  term  fluctuations.  The  layout  of  the  instrument  arrangement  is  shown  in 
Fig.5.6 
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Figure  5.4:  Velocity  Distribution  at  I=Idry 
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Figure  5.5:  Instrument  CeJibrations 
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Figure  5.6:  Instrument  Connections 

5.3  Experimental  Arrangement  and  Measurements 

The  selections  of  test  conditions  ideally  should  fulfill  the  following  criteria: 

1.  Inlet  and  beach  configurations  representative  of  nature. 

2.  Waves  and  current  conditions  in  conformance  with  the  basic  assumptions  in  the 
analytical  model. 

3.  Test  conditions  covering  a wide  spectrum. 

Due  to  the  limitations  of  the  basin  facility,  compromises  are  made.  The  beach  slope 
was  restricted  to  1:13.3  which  is  somewhat  steeper  than  most  of  the  natural  beaches  along 
the  east  coast  of  the  United  States.  The  inlet  dimensions  were  such  that  the  depth  to  width 
ratio  was  1:32,  which  was  common  for  the  inlet  along  Barrier  Island  and  the  model  size 
when  scaled  up  100  times  roughly  represents  a typical  natural  inlet.  The  test  range  on 
wave  conditions  was  unfortunately  rather  limited  owing  to  the  inability  of  the  wave  maker 
to  produce  clean  large  amplitude  waves.  Thus,  the  selected  incident  wave  amplitudes  were 
all  quite  small.  The  range  of  wave  period  (or  length)  variation,  was  also  limited  because  of 
the  size  of  the  basin  and  the  strong  reflection  due  to  the  beach  and  the  current  filter  at  the 
end  of  the  inlet. 

Six  different  cases  were  tested.  The  corresponding  conditions  are  given  in  Table  5.1. 
The  Ursell  number  is  an  indicator  of  wave  nonlinearity.  The  values  selected  here  were  to 
insure  that  the  Stokes  theory  remains  valid  at  least  in  the  constant  depth  region.  The  ratios 
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of  current  velocity  to  wave  celerity  at  the  exit  of  the  inlet  ranged  from  1:5.5  to  1:7.7  which 
were  within  the  range  that  the  basic  assumption  in  the  wave  energy  equation  is  not  violated. 

Measurements  were  performed  along  sections  parallel  to  the  y-axis  as  shown  in  Fig. 5. 3, 
waves  were  measured  along  six  sections  at  i = 4,  16,  26,  31,  35  and  37.  Current,  owing  to  the 
rapidly  diminishing  strength  from  the  inlet  exit,  were  measured  along  three  cross-sections 
only  at  i = 26,  31  and  35.  At  each  station,  relocation  at  a minimum  of  five  vertical  locations 
were  measured,  each  with  a time-averaged  value  of  3-5  minutes.  Since  the  inlet  discharge  is 
perpendicular  to  the  y coordinate,  the  velocity  component  in  the  y-direction,  V,  was  small. 
No  attempt  was  made  to  measure  this  component. 


Table  5.1:  Test  Condition 


Qoim^/s) 

Ho{m) 

T{sec) 

aoko 

koho 

Ur 

0.0097 

0.87 

0.027 

2.03 

0.25 

0.0056 

0.0088 

1.16 

0.015 

1.27 

0.58 

0.0097 

1.36 

0.013 

1.02 

0.99 

0.0123 

0.87 

0.039 

2.03 

0.32 

0.0076 

0.0095 

1.16 

0.017 

1.27 

0.62 

0.0097 

1.36 

0.013 

1.02 

0.99 

in  which 

Qo'-  Total  discharge  from  inlet; 

Hq:  Incident  wave  height  at  water  depth  ho; 

T : Wave  period; 

uq:  Incident  wave  amplitude  = Hq/2; 

ko:  Incident  wave  number  at  water  depth  ho; 

ho:  Maximum  water  depth  = 0.369  m.; 

TT  T 2 

Ur'.  Ursell  number  = ■ ° in  which  Lo  is  the  wave  length  at  depth  ho- 


5.4  Experimental  Results  and  Numerical  Comparisons 


The  experimental  results  were  presents  here  together  with  numerical  computations  to 
offer  one  to  one  comparisons.  Since  only  wave  amplitudes  and  x-component  of  the  mean 
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current  velocities  were  measured  in  the  experiments,  the  comparisons  were  also  restricted  to 
these  two  parameters.  Other  flow  properties  such  as  current  flelds,  wave  height  distributions 
and  stream  lines  were  presented  in  numerical  computations  only.  The  conditions  were  steady 
state,  which,  in  the  numerical  computation,  normally  took  500  iterations  with  a time  step  of 
0.5  second.  The  grid  system,  both  mesh  size  and  grid  layout,  was  identical  for  the  numerical 
and  physical  models. 

The  numerical  model  requires  that  the  mean  water  level  at  the  offshore  boundary  be 
specified.  Due  to  the  small  aoA:o  values  used  in  the  experiments,  the  actual  water  surface 
elevation  under  the  influence  of  set  down  is  difficult  to  established.  Therefore,  the  sensitivity 
of  the  numerical  model  due  to  slight  variations  of  mean  water  level  was  tested.  The  results 
of  current  velocity  and  wave  height  profile  along  the  onshore  direction  are  shown  in  Fig. 5. 7 
and  5.8,  for  rj  = 0 and  r}  = ~ 2sinh2fc£)  down  value  of  linear  wave  theory).  The  other 

input  conditions  are:  Hq  = I cm,  T ^ I sec.,  Qq  = 0.0019  m®/sec,  m = 13,  Ax  = 0.12 

® 0.18  m and  the  inlet  located  at  j = 3 and  4.  The  results  shown  in  these 

figures  are  for  those  on  the  plane  beach  ( j = 2)  and  along  the  inlet  ( j = 3).  The  differences 
due  to  the  two  boundary  conditions  are  so  small  that  in  the  subsequent  computations,  the 
offshore  boundary  water  level  displacement  is  set  equal  to  zero. 

The  grid  size  in  the  x-direction  controls  the  accuracy  and  the  rate  of  convergency.  A 
value  of  Ax  = 0.12  m was  selected  after  several  numerical  trails.  Figure  5.9  compares  the 
velocity  and  wave  profiles  for  the  same  testing  condition  with  Ax  = 0.12  m and  0.06  m, 
respectively.  For  Ax  = 0.06  m the  computational  time  doubles  that  for  Ax  = 0.12  m.  The 
maximum  difference  between  the  results  is  less  than  8 %. 

5.4.1  Wave  Shoaling  on  Plane  Beach 

The  first  case  compared  was  the  wave  shoaling  on  plane  beach.  The  test  conditions 
were.  Hq  = 1.23,  0.95,  0.97  cm  and  T — 0.87,  1.16,  1.36  sec.,  respectively.  The  results  are 
shown  in  Fig. 5. 10  with  calculations  from  both  linear  shoaling  and  the  full  equation  (2.138). 
Both  calculations  underestimated  the  shoaling,  particularly  for  longer  period  waves.  The 
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Figure  5.7:  Velocity  Variations  along  the  Beach,  — »7(l,y)  = 0;  o »7(l,i)  = -a^/{2  sinh  kh) 
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Figure  5.8:  Wave  Height  Variations  along  the  Beach, 

= -aV(2sinhA:h) 
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Figure  5.9:  Comparisons  between  m = 43  with  Ax  = 0.12  m, ; 2ind  m = 83  with  Ax 

= 0.06  m,  - - - at  II  = 15,30  for  m =43  and  12  = 29,59  for  m = 83 
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full  equation  actually  predicts  smaller  shoaling  values  than  the  linear  equation.  This  is  a 
main  shortcoming  of  wave  equation  presented  in  the  present  form  as  it  is  derived  by  using 
Stokes  parameter  to  rank  the  order  of  magnitude.  After  wave  breaking,  the  full  equation  is 
more  consistent  with  the  trend  of  the  experimental  data. 

5.4.2  Flow  Field  due  to  Inlet  on  Sloping  Beach  with  no  Waves 

The  second  set  of  experiments  was  conducted  with  inlet  flow  discharging  onto  a beach 
but  no  incoming  waves.  Two  cases  with  discharge  equal  to  0.0056  m^/sec  and  0.0076 
m^/sec  were  examined.  The  results  of  velocity  distributions  were  given  in  Fig. 5. 11  and 
5.12,  respectively.  The  computed  velocity  profiles  based  upon  the  numerical  model  were 
fitted  to  the  experimental  data  at  the  offshore  outer  cross  sections  i = 26  and  31  by  adjusting 
the  bottom  friction  coefficient  (/)  and  the  eddy  viscosity  coefficient  (cj.  The  reason  for 
selecting  the  outer  cross  sections  for  the  matching  is  that  both  frictional  and  diffusive  effects 
will  not  be  fully  recognized  until  the  flow  leaves  the  inertial  dominated  phase.  The  values 
for  / and  Cc  obtained  from  obtained  from  the  best  data  fit  are  as  follows: 

For  Qo  = 0.0056m^/sec 

/ = 0.005  and  e,.  - 0.003m^/sec 
For  Qo  = 0.0076m^/sec 

/ = 0.003  and  t,.  = 0.005m^/sec 

The  computed  velocity  distributions  are  in  good  agreement  with  the  data  both  in  magni- 
tude and  spreading  even  though  the  magnitude  is  somewhat  force  fitted.  Both  experimental 
data  and  calculated  profile  revealed  the  flow  reversal  characteristics  at  the  outer  edges  of 
the  jet  as  the  result  of  an  adverse  force  gradient  when  the  inlet  flow  enter  a sloping  beach. 

The  magnitudes  and  the  trend  of  variations  of  both  the  bottom  friction  coefficient  and 
the  viscosity  coeflScient  obtained  here  are  reasonable.  The  bottom  friction  coefficient  is 
smaller  for  the  larger  discharge  case.  This  is  because  the  Reynolds  numbers  in  the  inlet 
are  in  the  order  of  4,000  to  5,000  and  the  friction  coefficient  in  this  range  is  expected  to 
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decrease  as  the  velocity  increases.  The  eddy  viscosity  coefficient  is  a flow  property  affected 
by  the  mean  velocity  gradients;  one  expects  the  higher  the  velocity  gradient  the  larger  the 
eddy  viscosity.  Therefore,  the  values  obtained  here  are  consistent. 

5.4.3  Inlet  on  Sloping  Beach  with  Waves 

The  ultimate  purpose  of  the  present  study  is  the  development  and  verification  of  a 
numerical  model  that  is  capable  of  solving  the  kinematic  and  dynamic  flow  properties  of  an 
inlet-beach-wave  system.  Therefore,  the  measurements  of  this  category  are  of  paramount 
importance  to  test  the  validity  of  the  numerical  model. 

As  mentioned  before,  a total  of  six  different  combinations  of  current  and  wave  conditions 
were  tested.  The  effect  of  waves  on  current  field  modifications  and  influence  of  current  on 
wave  properties  are  discussed  here  separately. 

(1)  Current  Field  Modification 

The  experimental  results  of  current  profiles  for  the  six  cases  are  given  in  Fig.5.13  to 
5.24.  The  calculated  values  from  the  numerical  model  are  also  plotted.  The  solid  line  is 
the  computed  profile  with  the  same  input  conditions  as  the  experiment.  The  dotted  line  is 
the  corresponding  profile  without  input  waves.  In  these  computations,  the  eddy  viscosity 
coefficients  were  adopted  from  previous  experiments  where  no  waves  were  imposed  and  are 
assumed  to  be  a function  of  inlet  discharge  only.  The  rationale  behind  such  a choice  is 
based  on  the  physical  reasoning  that  the  lateral  mixing  caused  by  turbulence  generated  by 
current  gradient  will  not  be  significantly  affected  by  an  oscillatory  wave  field  and  on  the 
observed  fact  that  the  computed  velocity  profile  is  not  sensitive  to  small  changes  of  this 
coefficient.  The  bottom  friction  coefficient,  on  the  other  hand,  by  virtue  of  its  definition 
(see  Eq.(2.136))  is  expected  to  vary  with  the  wave  induced  bottom  velocity.  Therefore,  the 
same  method  employed  in  the  no  wave  case  is  used  here  by  matching  the  profiles  at  outer 
cross  sections  through  adjusting  the  bottom  friction  coefficient.  Here  the  values  of  / and 
€c  deduced  from  experimental  and  numerical  comparison  are  summarized  in  the  following 
table.  The  variations  and  characteristics  of  / will  be  discussed  later. 
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case 

Table  5.2:  The  values  of  and  / in 
discharge  (m^/sec)  Hq  (cm)  T (sec) 

the  Tests 
€c(m^/ sec) 

f 

a 

0.0056 

0.97 

0.87 

0.003 

0.011 

b 

0.0056 

0.88 

1.16 

0.003 

0.013 

c 

0.0056 

0.97 

1.36 

0.003 

0.015 

d 

0.0076 

1.23 

0.87 

0.005 

0.012 

e 

0.0076 

0.95 

1.16 

0.005 

0.010 

f 

0.0076 

0.97 

1.36 

0.005 

0.0105 

Overall  the  experimental  date  agreed  well  with  the  laboratory  results.  Again,  the 
agreement  on  the  current  magnitude  was  somewhat  artificial  because  of  the  adjustment 
made  on  the  friction  coefficients.  Therefore,  the  agreement  of  the  shape  of  the  distribution 
is  more  indicative  of  the  validity  of  the  numerical  model.  For  all  the  cases  tested,  the 
influence  of  waves  on  the  current  field  was  small.  In  general,  an  opposing  wave  field  reduces 
the  centerline  velocity  of  the  jet  flow  from  the  inlet  but  causes  the  jet  to  spread  wider. 
Return  flows  appear  at  the  outer  edges  of  the  jet.  In  the  far  field,  a weak  undulating 
current  field  outside  the  jet  zone  is  developed.  This  current  field,  undulating  in  the  lateral 
direction,  becomes  more  prominent  for  higher  discharge,  an  indication  that  it  might  be 
strongly  influenced  by  the  lateral  mixing.  To  reduce  the  possibility  of  numerical  instability 
on  the  unsteady  nature,  the  numerical  iterations  were  raised  up  to  900  with  practically  the 
same  results. 

As  stated  before  the  current  strength  at  the  center  portion  of  the  jet  is  reduced  under  the 
influence  of  the  opposing  waves.  When  the  wave  height  increases,  the  velocity  distribution 
can  no  longer  maintain  bell  shaped  and  the  center  portion  becomes  concaved  as  clearly 
shown  m cases  a and  d in  the  numerical  calculations.  Experimental  results  also  displayed 
such  behavior,  only  more  pronounced.  In  case  /,  for  instance,  the  experimental  results 
showed  a very  prominent  concave  shape  whereas  the  numerical  prediction  only  revealed  a 
very  mild  indentation. 

Another  experimental  observation  is  the  seemingly  strong  three  dimensional  flow  pat- 
tern along  the  edges  of  the  jet  where  return  flow  developed  near  the  surface  layer. 

To  illustrate  the  change  of  the  complete  velocity  field,  contour  plots  of  current 
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magnitude  for  all  the  cases  are  shown  in  Figs. 5.25  to  5.34.  All  the  contour  lines  are  expressed 
in  non-  dimensional  velocity  defined  as  u/{Qq/A)  where  u = is  the  magmitude 

of  velocity,  Qq  is  the  inlet  discharge  and  A is  the  inlet  cross-sectional  area.  Figure  5.25 
shows  a case  of  an  inlet  of  Qq  = 0.0056m^/sec  discharging  into  a basin  of  flat  bottom  with 
constant  water  depth  of  0.045  m and  bottom  friction  coffficient  equal  to  zero.  In  this  case, 
the  flow  behaves  more  or  less  like  a two-dimensional  submerged  jet.  The  core  of  the  jet 
contains  a clear  out-shot  tongue.  The  extent  of  the  tongue  is  an  indicator  of  a flow  regime 
dominated  by  inertial  force.  The  core  ends,  of  course,  where  the  zone  of  diffusion  reaches 
the  centerline  of  the  jet.  Downstream  from  this  section,  the  magnitude  of  the  velocity  di- 
nunishes  as  the  jet  expands  laterally  due  to  entrainment.  The  velocity  distribution  is  bell 
shaped  with  discernable  similarity  profiles  as  the  jet  continually  expands  into  downstream 
section.  Figure  5.26  shows  a similar  case  with  the  exception  that  bottom  friction  is  present 
(/  = 0.01).  The  presense  of  bottom  friction  greatly  retards  the  jet  flow  and  results  in  much 
rapid  lateral  mixing.  The  velocity  vectors  are  also  plotted  for  both  cases. 

For  the  situation  that  a jet  expands  into  a sloping  beach,  such  as  that  in  the  present 
experiment,  it  appears  that  the  jet  contracts  first  then  expands  gradually  with  pulsatory 
nature.  It  is  apparent  that  similarity  velocity  profile  can  not  be  maintained.  The  pulsatory 
nature  of  the  jet  profile  could  be  the  consequence  of  two  competing  mechanisms:  the  lat- 
eral diffusion  mechanism  prompts  the  jet  to  expand  but  the  cross-sectional  expansion  due 
to  water  depth  deepening,  on  the  other  hand,  creates  vertical  vortex  and  causes  the  jet 
expansion  in  the  vertical  direction  at  the  expense  of  lateral  contraction. 

The  apparent  effects  of  opposing  waves  are  the  retardation  of  the  advancing  jet  and  the 
wider  spreading  of  the  current  field.  The  over  all  effects  on  the  jet  appear  to  be  mild. 

As  mentioned  earlier,  the  values  of  bottom  friction  coefficient  were  obtained  by  the  best 
data  fitting  at  outer  cross  sections  i = 26  and  31.  These  values  are  list  in  Table  5.3  together 
with  other  pertinent  parameters.  Since  in  the  present  model  the  bottom  shear  stress  is 
computed  by  a pair  of  simplified  quadratic  equations  (Eqs.(2.118)  and  (2.119))  in  which 
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the  wave  obital  velocity  appears  only  as  a linear  quantity,  it  is  expected  that  the  associated 
friction  coefficient  should  be  affected  by  the  relative  magnitude  of  this  wave  orbital  velocity 
to  the  ambient  current  velocity.  The  values  of  / are  plotted  against  u^jUi  in  Fig.5.35  in 
which  Ui  is  the  discharge  velocity  at  inlet  and  the  definition  of  is 

o,u; 

Um  = 

sinh  kh 

where  a,-  is  the  wave  amplitude  at  inlet  given  by  linear  shoaling  and  h is  the  still  water  depth 
at  inlet,  k is  wave  number  at  inlet  with  no  current.  It  appears  that  / in  the  simplified  bottom 
stress  equation  is  linearly  proportional  to  The  results  so  obtained  are  completely 

consistent  with  the  physical  phenomenon. 


Table  5.3:  Bottom  Friction  Coefficients  Used  in  the  Model 


case 

T{sec) 

oo  (cm) 

Um  (cm/s) 

Uq  (cm/s) 

^yti/Uq 

U 

a 

0.87 

0.485 

6.58 

8.64 

0.76 

0.011 

b 

1.16 

0.440 

6.20 

8.64 

0.72 

0.013 

c 

1.36 

0.485 

6.92 

8.64 

0.80 

0.015 

d 

0.87 

0.615 

8.35 

11.73 

0.71 

0.012 

e 

1.16 

0.475 

6.69 

11.73 

0.57 

0.010 

f 

1.36 

0.485 

6.92 

11.73 

0.59 

0.0105 

(2)  Wave  Field  Modification 

In  contrast  to  the  good  agreement  between  measured  and  computed  velocity  field,  the 
wave  amplitudes  were  mostly  underestimated  by  the  numerical  model.  This,  as  explained 
earlier,  is  mainly  due  to  the  limitation  of  the  wave  equation  in  the  form  used  in  the  present 
model,  i.e.,  to  include  the  effects  of  refraction,  diffraction  and  shoaling  on  the  same  order 
of  magnitude. 

In  the  experiment,  owing  to  the  finite  basin  size,  the  presence  of  inlet  and  current 
influenced  the  wave  field  of  the  entire  basin  including  the  boundary  wave  heights.  In  the 
numerical  computation,  the  measured  wave  amplitudes  at  the  offshore  boundary  are  used 
as  input.  Table  5.4  provides  the  boundary  heights  as  measured.  The  centerline  of  the 
basin  is  in  between  grid  20  and  21.  For  normal  incident  waves,  the  field  is  symmetrical. 
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Table  5.4:  Incident  Wave  Height 


case 

Qq  (m^/sec) 

T (sec) 

Hq  (cm) 

Hn 

^18 

^19 

H20 

a 

0.0056 

0.87 

0.97 

1.00 

1.03 

1.06 

1.07 

b 

0.0056 

1.16 

0.88 

0.88 

0.89 

0.92 

0.965 

c 

0.0056 

1.36 

0.97 

0.97 

0.98 

0.99 

0.99 

d 

0.0076 

0.87 

1.23 

1.27 

1.33 

1.42 

1.44 

e 

0.0076 

1.16 

0.95 

0.97 

0.98 

1.05 

1.08 

f 

0.0076 

1.36 

0.97 

0.99 

0.99 

1.13 

1.15 

Therefore, 

we  have  H21  = 

H20,  H22  - 

Hig  and  so  on. 

The  wave 

heights 

of  unspecified  grid 

are  considered  to  be  equal  to  Hq,  where  Hq  is  the  original  incident  wave  height. 


Figures  5.36-5.41  give  the  measured  data  plotted  against  the  numerical  values.  Both 
numerical  and  experimental  results  revealed  that  the  wave-current  interaction  has  a very 
prominent  effect  on  the  local  wave  height,  particular  in  the  zone  directly  influenced  by 
the  jet.  The  counter  current  in  the  center  of  the  jet  obviously  enhances  the  wave  height 
whereas  the  return  current  near  the  edge  of  the  jet  causes  the  waves  to  diminish  somewhat. 
Therefore,  the  shape  of  the  wave  height  profile  is  very  similar  to  that  of  the  current.  Figure 
5.42  shows  the  variations  of  the  non-dimensional  wave  heights  along  the  centerline  of  the 
basin  for  the  six  cases  tested.  As  can  be  seen,  the  combined  shoaling  and  current-wave 
interaction  causes  the  waves  to  increase  their  amplitudes  significantly  as  they  progress 
toward  the  inlet.  The  numerical  results  showed  more  than  doubling  the  amplitude  for  cases 
of  strong  current.  The  experimental  results  were  even  higher,  with  amplification  reaching 
as  high  as  threefold. 


The  agreement  between  experimental  and  numerical  results  is  considered  fair  for  weak 
current  and  in  deep  water  and  becomes  progressively  worse  with  longer  incoming  wave 
period,  shallower  depth  and  stronger  counter  current.  The  breaking  criterion  used  in  the 
numerical  model  proves  to  be  inadequate  under  strong  current  condition.  Experimental 
evidence  showed  that  waves  break  earlier  than  predicted. 
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Figure  5.10:  Wave  Height  Variations  along  the  Beach  (no  current), Wave  Eq.  (2.138); 

Linear  Shoaling;  x Measured  Data 
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Figure  5.11:  Velocity  profiles  for  Qo  = 0.0056m^/sec.  (no  waves) 
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Figure  5.12:  Velocity  profiles  for  Qo  — 0.0076m^/sec.  (no  waves) 
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Figure  5.13:  Variation  of  x-Component  Velocity.  Case  a:  Qq  = 0.0056m^/«;  T = 0.87  sec. 
Hq  — 0.97  cm. Numerical  Model;  o Measured  Data 
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Figure  5.14:  Comparisons  of  Calculating  Velocities  between  with  and  without  Waves.  Case 
a:  Qo  = 0.0056m*/s,  T = 0.87  sec.;  Hq  = 0.97  cm.  with  Waves; no  Waves 
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Figure  5.15;  Variation  of  x-Component  Velocity.  Case  b:  Qo  = 0,0056m^/s;  T = 1.16  sec. 
Hq  = G.88  cm.  Numerical  Model;  o Measured  Data 
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Figure  5.16:  Comparisons  of  Calculating  Velocities  between  with  and  without  Waves, 
b:  Qo  = 0.0056m3/s,  T = 1.16  sec.;  Hq  = 0.88  cm.  with  Waves;  - - - no  Waves 


Case 
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Figure  5.17:  Variation  of  x-Component  Velocity.  Case  c:  Qq  = 0.0056m®/s;  T = 1.36  sec. 
He  = 0.97  cm.  Numerical  Model;  o Measured  Data 
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Figure 

c:  Qo  - 
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5.18:  Comparisons  of  Calculating  Velocities  between  with  and  without  Waves.  Case 
C.0056m^/s,  T = 1.36  sec.;  Hq  = 0.97  cm.  with  Waves; no  Waves 
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Figure  5.19.  Variation  of  x-Component  Velocity  Case  d:  Qq  = 0.0G76m'’'/s;  T = 0,87  sec.; 
Hq  = 1.23  cm.  NumericaJ  Model;  o Measured  Data 
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Figure  5.20:  Compaxisons  of  Calculating  Velocities  between  with  and  without  Waves.  Case 
d:  Qo  = 0.0076m’/fi,  T = 0.87  sec.;  Ho  = 1.23  cm.  with  Waves; no  Waves 


Ill 


o 

o 


o 


o 


Figure  5.21:  Variation  of  x-Component  Velocity.  Ceise  e:  Qq  = 0.0076m^/s;  T = 1.16  sec. 
Hq  = 0.95  cm. Numerical  Model;  o Measured  Data 


U (CM/SEC)  U (CM/SEC)  U (CM/SEC) 
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Figure  5.22:  Comparisons  of  Calculating  Velocities  between  with  and  without  waves.  Case 
G-  Qo  — 0.0076m^/s,  T = 1.16  sec.;  Hq  = 0.95  cm.  with  Waves; no  Waves 
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Figure  5.23:  Variation  of  x-Component  Velocity.  Case  f:  Qq  — 0. 0076m® /s;  T 
^0  — 0.97  cm.  Numerical  Model;  o Measured  Data 


= 1,36  sec.; 
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Figure  5.24:  Comparisons  of  Calculating  Velocities  between  with  and  without  Waves.  Case 
f-  Qo  — 0.0076m^/«,  T = 1.36  sec.;  Ho  = 0.97  cm. with  Waves;  - - - no  Waves 
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Figure  5.25:  Velocity  Field  and  Contour  Lines  in  Flat  Bottom. 
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Figure  5.26:  Velocity  Field  and  Contour  Lines  in  Flat  Bottom,  f — 0.01 
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Figure  5.27:  Contour  Plot  of  Velocity.  Case  :.  Qo  — 0.0056m^/s  (no  waves) 
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Figure  5.28:  Contour  Plot  of  Velocity.  Ceise  a:  Qo  = 0.0056m®/s;  T = 0.87;  Hq  = 0.97  cm 
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Figure  5.29:  Contour  Plot  of  Velocity.  Case  b:  Qo  — 0. 0056m® /s;  T — 1.16,  Hq  0.88  cm 


120 


Figure  5.30;  Contour  Plot  of  Velocity.  Caise  c:  Qo  = 0.0G56m^/s;  T = 1.36;  Ho  = 0.97  cm 
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Figure  5.31:  Contour  Plot  of  Velocity.  Case  ii:  Qo  = 0.0076m’/s  (no  waves) 
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Figure  5.32:  Contour  Plot  of  Velocity.  Case  d:  Qo  = O.OOTGm^/s;  T = 0.87;  Hq  = 1.23  cm 
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Figure  5.33:  Ccptour  Plot  of  Velocity.  Case  e:  Qo  = 0.0076m^/s;  T — 1 16;  ifo  ” ^-95  cm 
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Figure  5.34:  Contour  Plot  of  Velocity.  Caise  f:  Qo  = O.CX)76m*/s;  T — 1.36;  Ho  = 0.97  cm 
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Figure  5.35:  Bottom  Friction  Coefficients  / as  Function  of  u^/U. 
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Figure  5.36:  Wave  Amplitude  Variations.  Case  a:  Qo  = 0.0056m^/s;  T = 0.87  sec.;  Hq  =- 
0.97  cm.  Numerical  Model;  A Data 
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Figure  5.37:  Wave  Amplitude  Variations.  Case  b:  Qq  = 0.0056m^/s;  T = 1.16  sec.,  Hq  = 
0.88  cm. Numerical  Model;  A Data 
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Figure  5.38;  Weve  Amplitude  Variations.  Case  c:  Qq  = 0.0056m^/s;  T = 1.36  sec.;  Hq  = 
0.97  cm.  Numerical  Model;  A Data 
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Figure  5.39;  Wave  Amplitude  Variations.  Ceise  d:  Qo  — 0.0076m^/s.  T — 0.87  sec.;  Ifo  = 
1.23  cm.  Numerical  Model;  A Data 
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Figure  5.40:  Wave  A.mplitude  Variations.  Case  e:  Qq  — O.OOTGnr/s;  T = 1.16  sec.;  Hq  = 
G.05  cm. Numeiical  Model;  A Data 
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Figure  5.41:  Wave  Amplitude  Variations.  Case  f:  Qq  0,0076fn^/s;  T --  1.38  eec.;  Ho  ~ 
0.97  cm.  Numerical  Model;  A Data 
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Figure  5.42:  Wave  Amplitude  Variations  along  Centerline 


CHAPTER  6 

CONCLUSIONS  AND  RECOMMENDATIONS 


The  main  objective  of  the  present  study  is  to  develop  a numerical  model  capable  of 
predicting  the  hydrodynamics  in  the  vicinity  of  an  ocean  inlet,  which  is  a complicated 
system  owing  to  the  interactions  of  spatially  non-uniform  current  with  waves  over  changing 
topographies.  The  study  entails  three  major  endeavors:  theoretical  formulation,  numerical 
modeling  and  experimental  verification. 

The  theoretical  formulation  includes  the  derivations  of  four  basic  equations:  the  depth- 
averaged  continuity  equation,  two  depth-averaged  momentum  equations  and  the  wave  equa- 
tion. The  derivation  of  the  first  three  equations  is  straight  forward  and  follows  closely  the 
procedures  employed  by  Ebersole  and  Dalrymple  (1979).  The  wave  equation,  on  the  other 
hand,  is  extremely  complex  and  has  to  be  dealt  with  care  and  tailored  to  the  physical  system 
it  intended  to  represent. 

The  wave  equation  derived  here  retained  the  important  effects  of  refraction,  diffraction, 
current-wave  interaction  and  energy  dissipation.  The  wave  refiection  effect  is  not  considered 
under  the  assumption  of  slowly  varying  topography.  Two  different  method  - variational 
principle  and  perturbation  method  - are  used  which  yield  identical  final  form  of  the  wave 
equation.  One  primary  magnitude  scale,  e,  and  two  auxiliary  magnitude  scales,  /z  and  8 are 
introduced,  e is  the  commonly  known  Stokes  wave  steepness  parameter;  /z  is  a spatial  scale 
parameter  representing  the  relative  horizontal  spatial  gradient  over  one  wave  length  and 
finally,  5 is  a strength  parameter  indicating  the  relative  strength  of  mean  current  to  wave 
celerity.  In  the  present  model,  the  wave  equation  is  subject  to  the  following  constraints: 
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which  imply  slowly  varying  topography,  strong  mean  current  and  small  incident  wave  angle. 
The  final  form  of  the  equation  is  consistent  to  O(e^)  and  contains  all  pertinent  terms  to  that 
order.  Therefore,  consistency  (or  homogeneity)  and  completeness  are  the  emphasis  in  the 
present  derivation. 

The  equations  form  a set  of  nonlinear  partial  differential  equations  dependent  upon 
time  and  space.  The  numerical  method  developed  to  solve  these  equations  is  based  upon 
Crank-Nicolson  finite  difference  scheme  in  a double  spatial  sweep  form.  A two  time  level 
implicit  model  developed  by  Sheng  (1981)  is  incorporated  into  both  the  x-sweep  and  the 
y-sweep.  To  implement  the  model  requires  specification  of  the  bottom  friction  coefficient, 
lateral  diffusion  coefficient  as  well  as  the  boundary  conditions  of  waves  and  current. 

Performance  of  the  model  is  evaluated  against  existing  numerical  models  including  cases 
of  wave  set-up  on  a plane  beach,  longshore  current  distributions  due  to  an  oblique  waves 
on  a plane  beach,  wave-induced  nearshore  circulation  over  an  irregular  topography  and 
wave  diffraction  over  a circular  shoal.  Examples  are  also  given  to  demonstrate  the  model’s 
capability  of  providing  nearshore  current  and  wave  interaction  for  an  inlet-beach  system, 
which  is  the  ultimate  goal  of  this  study. 

Finally,  laboratory  experiments  were  conducted  to  calibrate  and  verify  the  numerical 
model  and  to  evaluate  the  two  empirical  coefficients  - bottom  friction  and  lateral  diffusion 
required  as  input  to  the  model.  The  laboratory  conditions  are  intended  to  represent 
the  typical  coastal  inlets  along  the  eastern  seaboard  of  the  United  States.  Owing  to  the 
limitation  of  the  facility,  only  a limited  range  of  inlet  current  strength  and  wave  conditions 
were  tested.  By  adjusting  diffusion  and  friction  coefficients,  the  numerical  model  could 
duplicate  the  current  profiles  satisfactorily.  However,  the  model  invariably  underpredicts  the 
wave  shoaling  against  current.  This  is  partially  due  to  the  inherent  weakness  of  the  governing 
wave  equation  and  partially  due  to  the  inadequate  breaking  criterion.  The  lateral  diffusion 
coefficient  and  the  bottom  friction  coefficient  so  evaluated  by  force  fitting  the  current  profiles 
are  found  to  be  reasonable  and  their  variations  can  be  explained  on  plausible  physical 
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reasoning.  Both  physical  and  numerical  models  reveal  that,  in  the  range  of  commonly 
encountered  inlet  conditions,  the  effects  of  current  on  waves  are  much  more  pronounced 
than  the  other  way  around.  In  the  physical  model,  wave  amplification  as  high  as  threefold  is 
observed  while  numerical  prediction  also  yields  values  more  than  doubling  the  incident  wave 
heights.  Topography  is  still  the  dominant  factor  governing  the  current  pattern.  The  effects 
of  bottom  friction  appear  to  be  more  pronounced  than  the  lateral  diffusion  in  retarding  the 
inlet  flow  and  causing  wide  spreading  of  the  jet. 

The  numerical  model  developed  here  should  prove  to  be  a useful  tool  for  better  un- 
derstanding of  inlet  dynamics  and  as  a basis  to  develop  a sediment  transport  model.  The 
emphasis  on  the  inclusion  of  wave  diffraction  effects  necessitates  compromising  on  wave 
shoaling  prediction  as  well  as  imposing  stricter  restrictions  on  the  incoming  wave  angle. 
Although  the  limit  on  mild  slope  criterion  cannot  be  clearly  established,  it  is  believed  that 
the  model  is  able  to  accommodate  common  prototype  nearshore  topographies  without  nu- 
merical errors. 

Although  the  numerical  model  seems  to  yield  sensible  results,  the  verification  is  inade- 
quate and  incomplete  for  two  reasons:  most  of  existing  numerical  models  used  as  bench  mark 
here  are  themselves  not  adequately  verified  and  the  physical  model  experiments  conducted 
in  this  study  are  too  restricted  due  to  facility  and  instrument  limitations. 

It  is  recommended  here  that  the  first  and  foremost  teisk  is  to  expand  the  laboratory 
experiments  to  verify  the  numerical  model  with  particular  emphasis  on  current  measure- 
ments. Field  data  on  current  in  the  vicinity  of  inlet  should  also  prove  to  be  invaluable  as 
we  now  have  the  preliminary  tool  to  interpret  them. 


APPENDIX  A 

WAVE  EQUATION  DERIVED  BY  USING  PERTURBATION  METHOD 


Appendix  A gives  a derivation  of  the  wave  equation  based  on  perturbation  method.  It 
follows  Chu  and  Mei  (1970),  but  under  the  consideration  of  strong  current.  On  the  other 
hand,  instead  of  unsteady  time  scale,  steady  time  scale  is  considered,  it  means  that  no  slow 
derivative  is  related  to  time. 

Governing  equations  read 


V^<j>  = 0 —ho  < z < T) 

<f>z  = Vt  + ^ h<l>  • ^ hV  z = rf 

+ I I V<^  P=  0 Z = T) 


4>z  + ^ h<f>  • V/,ho  = 0 z — —h 


where 


4>  = <f>c  + 4> , 
ri  = riQ  + r]' 

and  <j>c  and  r}o  represent  current  alone,  <t>'  and  r?'  are  wave  quantities  which  have  the  order 
— 0(^)>  ^0  is  still  water  depth.  V and  V;,  axe  again  the  three  dimensional  and  horizontal 
gradient  operators,  respectively.  As  in  chapter  2,  introduce  a shifting  coordinate  about  z 
axis,  which  define  as 

z'  - Z - T]o 

z -j-  Hq  = z'  h' 

A - A 

dz  dz' 
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Dropping  the  primes  on  z and  h for  convience,  governing  equations  are  replaced  by 


< 

to 

II 

o 

—ho  < Z < T]' 

(A.1) 

= v't  + • V/ifj'  + Vh<f>  • V/if/0 

z = r,' 

(A.2) 

<f>t  + givo  + f?')  + 5 1 p +f  1 V(/)  |2=  0 

z = rj' 

(A.3) 

4>z  + ^ h<t>  • ^ hh  = 0 

z — —h 

(A.4) 

in  which 


I |2  + I |2=|  p 

I \l=\  V^'  p +2V^'  • V<f>, 

^h4>c{^>y)  and  r)o{x,y)  are  0(1)  under  strong  current  assumption  as  did  in  chapter  2. 

= Ui  + Vj 

From  the  kinetic  free  surface  boundary  condition  z'  = z — rjo  we  get 


u • V/i  rjo  = 0 from  (A. 2) 

Collecting  0(1)  terms  in  (A.3)  gives 

yo  = (A.5) 

which  is  same  as  (2.64)  in  chapter  2;  finally  the  governing  equations  for  wave  motion  are 
given  by 

= 0 -ho<z<r)'  (A.6) 

9<pz  + 4>tt  + + l'V(f>-'V)  \v^\l=0  Z = r}'  (A.7) 

4>t  + gr)'  +l\V<f>\l^0  z^ri'  (A.8) 

+ V;,^  • V/ih  = 0 z=-h  (A.9) 

in  which  we  have  substituted  (A.3)  into  (A.2)  to  make  (A.7)  only  contains  function  <f>. 
Introducing  a slowly  varying  coordinates,  i.e. 


x = x + ^x  = X + X2 
y = + (iy^Yi  + Y2 
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where  fi  is  slowly  changing  parameter,  as  defined  in  chapter  2,  so  that 

d 2 d 

dx~ 

d _ d 2 d 

dy~^Wi'^^  W2 

Assuming  that  <f>  and  rj  have  the  forms 

3 n 

4>  = €-^<I>,{X2,Y2)  + Y.  E (A.IO) 

n=:l  m—~n 


3 n 

»?  = »?o(A2,y2)  + E E eV”’”‘Hn,V2,A2,z)e‘"‘^  (A.ll) 

n=l  m=—n 

where  0 is  phase  function  which  is  defined  under  assumption  that  waves  mainly  propagate  in 
and  let  complex  amplitude  absorbs  the  different  between  real  phase  and  imaginal 
main  phase.  In  supercripts  n is  the  orders  of  magnitudes  and  m is  the  harmonics,  m can 
not  be  great  than  n. 


.-f 


k\dx  — u)t 


ki  = k cos  6 
k2  = k sin  6 


We  have  assumed  that  incident  wave  angle  is  small,  therefore 

k^  = kl  + kl  = kl  + O(e^) 

Substituting  (A.IO)  and  (A.ll)  into  (A.6)  through  (A.9),  after  taking  Tailor  series  about 
2 = 0,  we  orginize  equations  as 


(^  - m^A:J)(/>("-”‘)  = ir(n.m) 

— h < Z < 0 

(A.12) 

^(n.m)  _ Q{n,m) 

z = —h 

(A.13) 

{g-^  - 

z = 0 

(A.14) 

2 = 0 

(A.15) 
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The  and  are  given  as  follows.  It  is  noted  of  that  | V h<t>c  ? is 

removed  from  | V<j>  |J. 

(a) 

Define  = [ | 

n m 


where  5 is  defined  as 


5 = 


1 ni  — ki  = ri2  — k2]  mi  = m2 

2 otherwise 


where  (nx  - ki)  + (ri2  - ^2)  = n and  mi  + m2  = m in  which  k{  - 0, 1, 2. 

The  common  term  e"  in  /(”■”*)  has  been  dropped.  Writing  in  details,  we  have 


/(o.o)  ^ 0 

/(1,0)  ^ 0 

/P.‘)  = 4*:f^PPV<‘“"’  + 2ri'P)4‘“‘>  + 2^P'“V1''‘’  + 24^’‘V1‘'°’ 
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Next  define 


Expanded  as 


Q{n,m) 


^(ni-2,mi) y{n2.m2)  ^ [ . . . j{n-3,m) 


Term  by  term,  yield 

q(°-°)  = 0 

= 0 

g(2.o)  ^ l<^(b-i)/(i.i)  + i0(bi)/(i.-i)  + 1 ^(i.i);(i,-i)  + ^(1,-1)  ^(1.1)  ] 

g(2.i)  ^ -.•c^/(2-^)  + i.-McX2/^'’'^  + ^^cr2/^^'’  + 

g(2.2)  = -2,-a;/(2.2)  + .-fci^,;,^/(2.2)  _ \k\cf>^hl)  f{l.l)  + 1<^(1.1)^(1.1) 

Q(«)  = -,„/(3,.)  + ^ 

+ ^^cr,4;’'>  +*;|(»<‘’-'l/(3'3>  + ^(3.3)/(l.-3)]  + i|4>.»/(3.«) 

2 

+^(l.0)/(2.1)  + ^(l.-l)/(2.2)  + ^(2,2)^(1.-1)  + ^(2.0) ^(1.1)  ] 

Finally,  expanding  Q(”-”‘)  in  Tailor  series  about  z = 0,  we  obtain 

Gr(n,m)  ^ _g(n,m)  _ j ^ _ 1 j ^ ^ ^ _ 

Writing  in  detail,  we  have 

+9<f.  + « I + }.>,  + (?  |(;>-3'"->)  + . . . 
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where 


^(n.m) 


2 ' 


fj(n-2,m2) 


The  expansions  of  are  as  follows: 


^(2.1)  ^ ^(1.0)^(1,1) 

^(2,-1)  ^ ^(1 -1)^(1, 0) 
f(2.2)  = 

2 

The  formations  of  are  given  as  follows: 

G!(o.o)  Q 

G(i.o)  = 0 

G(^.i)  = _q(i.i) 

Gft»)  = -Q(2,0)-,C.'»|s,^,  + <3 ](■.») -,(i.-i)|-a,V  + j^.  + QJ(‘.') 

G(=.»  = - ,C.0)(  + ,^,  + <3  J(^.0  - ,(M)|  + Q |(J,0)  _ ,(l,-i) 

I -4«V  + 9^1+13  ll*'’’  - -"V  + 9#«  + <3 

-wV  + 9*.  + <3 

(b) 

= -[  + <^cY2  - hY2<t>^Y^~^'"'^ 
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Thus 


5(1,0)  ^ 5(1,1)  _ 5(1,-!)  _ Q 

5(2,1)  ^ 5(2,2)  ^ 0 

5(3,1)  _ —tmkihx24>^^’^^  for  slip  bottom  boundary  condition 
or 

5(3,1)  — _|_  ,-j^^^y^2^(i,i)  ^ y^^(i,i)  _ non-slip  bottom 

boundary  condition 


The  non-slip  bottom  boundary  condition  is  given  in  Appendix  B after  following  Liu  (1986). 
(c)  f (”'"*) 


Writing  in  detail,  yield 


5(1,1)  _ 5(1,-!)  _ Q 

f(2,o)  -V;,Vc 

5(2,1)  ^ p{2,2)  ^ Q 

5(3,1)  _ ^ 2tki<f>x2  + 4>YiYi 

(d) 


/f(«.»«)  = + 1 I |2  ]('».”») 

9 2 

Expanding  above  equation  in  Tailor  series  about  z = 0,  we  obtain 


-[imuj4>-  i/jK-bma)  ^ i^(2,mi) 

9 2 g 2 g 

[im20J<l>  - If  ](”2-2.’"2)  + _ l/](«2-2,m2)  ^ . 

2 g 2 


Therefore 
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^(1.0)  ^ 0 

9 ^ 

= -i[,-o;,^+i/](i.-i) 

= _^;(!.0)  _ 1,(M)|,„^+  1;|(.,-.)  + 1,(1, _ ‘/](l.« 

- l/l^" 

9^9  ^ 

F(2.2)  = l[2,-o;^-  i/](2.2)  + 

9^9  2 

The  solutions  will  be  given  in  the  form  of  combinations  for  different  orders  and  har- 
monics. 


(1)  n = 1,  m = 0 


The  governing  equations  are  given  as 

=0  -h<  z <0 

= 0 z = -h 

= G(1-0)  = 0 z = 0 

^(1,0)  ^ ^(1.0)  ^0  2 = 0 


The  solutions  are  simply  obtained  as 


,,(1-°)  = 0 
41,0)  = ^(.,0)  ^ 0 


(2)  n = 1,  m = 1 


(A.16) 

(A.17) 


4>^zz^^  - A:f^(bi)  =0  -h  < z <0 

=0  z = -h 

94>i  ’ ^ 

94>^z'^^  - (w  - = 0 2 = 0 

_ ki^cXi  ] z = 0 


144 


Those  give 


^(1.1)  _ .•  9 _coshki{h  + z)  , 

(^v  . ; _ — a ^ — /g»v  ^ c.c 

l<Tl  coshKi/i 

+ c.c 

£i 


with  the  definition  of  the  apparent  dispersion  relation 


(A.18) 

(A.19) 


w = \/gk\  tanh  fcih  + ki<j>cX2 
= (Tl  + ki4>cX2 

in  which  a is  a complex  wave  amplitude  and  c.c  is  the  complex  conjugate;  ki  and  k2  are 
wave  numbers  in  x and  y direction,  respectively.  The  linear  dispersion  relation  discussed  in 
chapter  2 gives 


w = \/gk  tanh  kh  + ki<f>cX2  + k2<l>cY2 
= (T  + ki<f>cX2  + ^2<f>cY2 


Therefore,  we  have 


O’!  — O'  + k24>cY2  = o + 0(e) 


By  introducing  the  error  for  less  than  0(c*),  equation  (A.8)  can  be  revised  as 

^(M)  = 

2<7i  cosh  kh 


(A.20) 


It  is  important  to  note  that  when  we  insert  (A.20)  back  into  Laplace  equation,  we  get 

= {k^  - ~ o(e®)  (A.21) 


It  entirely  satisfies  governing  equation  at  order  of  0(e),  but  the  left  hand  side  will  be  induced 
at  the  order  of  O(e^). 

(3)  n = 1,  m = -1 


. g cosh  k(h  + z) 

t a -e  ^ + cc 

2<ti  cosh  kh 


+ c.c 


(A.22) 

(A.23) 


(4)  n = 2,  m = 0 


145 


- -^>0X2X2  - <i>cY2Y2 

-h  < z <0 

(^(2.0)  = 5(2.0)  _ -(j)^X2hx2  ~ <f>cY2^Y2 

z = —h 

^(2.0)  ^ ^(2,0)  ^ 0 

z = 0 

»?(2,0)  = i7(2.o)  = -^/(2.0)  + f 

z = 0 

The  following  solvability  condition  must  be  satisfied  by  first  above  three  equations: 

r F^^’°Uz  + 5(2-0)  ^ 1^(2, 0) 

J-h  g 

which  gives 

Vh{hV)  = 0 

(A.24) 

or 

hX2^X2  + hY2VY2  = 0 

(A.25) 

where  U = Vfi<f>c 

It  is  continuity  equation  for  order  of  0(1).  The  last  equation 

gives  r/(^’0)  as 

^(2.0)  _ 1 1 ^ 1_  0(e®) 

' 2sinh2A:/i  ( ^ 

(A.26) 

This  is  the  mean  water  set-down  related  to  the  shifted  coorinate  due  to  radiation  stresses. 

(5)  n = 2,  m = 1 

^(2-1)  - Jfc2^(2.i)  = 5(2,1)  ^ 0 

-h<  z<0 

4^-^)  = 5(2.1)  ^ 0 

z = —h 

- W242.I)  ^ (^(2,1) 

or  - «t2^(2.i)  = g'(2-^)  = 2t'cri<^cY2^Yi^^ 

z = 0 

^(2,1)  ^ ^(2,1)  ^ !^^(2.1)  + j^,Y2</>Y^’^^ 

2 = 0 

in  which  does  not  include  (^(^4) 
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The  possible  solutions  from  above  equations  are  and  equal  zeros.  However, 
»?(2,1)  = 0 contradicts  with  that  derived  in  chapter  2 which  t7(2,  1)  is  non-zero.  But, 
rj[2, 1)  will  not  be  involved  in  later  calculation  in  this  derivation,  it  does  not  affect  the  final 
formulation  of  wave  equation. 

(6)  n = 2,  m = 2 


- 4*2 = 0 

—h  < z < 0 

= 0 

z = —h 

- 4o;20(2.2)  _ C!(2,2)  qj. 

P0P)- 4^20(2.2) 

z = 0 

,,(2,2)  ^ ^(2.2) 

z = 0 

in  which 


i{uj  + tanh^  ] -|-  ikl<j>cX2<f>^^’^^^  ~ 


These  give  results  as 

exactly  the  same  as  the  second  order  Stokes  wave  results. 
(7)  n = 3,  m = 1 


(A.27) 

(A.28) 


— h < z < 0 

0(^-^)  = 5(3.1) 

z = —h 

^0(3.1)  — oi;^0(3.1)  = G(3.1) 

or 

<70^'^  - cr20(3.i)  = 

z = 0 
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where 


I - I 

-■/'“■“l  I + «&■'  + 1 - £<'■»)  I ■ ■ ■ 1 - f . . . I 


zero. 


The  summations  of  last  two  brackets  in  G!»{3.1) 

are  identically  equal 

The  last  term  in  right  hand  side  of  first  equation  is  due  to  — k\  = k\  = O(e^)  as 
mentioned  in  case  (2).  and  are  no  term  in  and  g(3,i)^  respectively. 

Again,  substituting  F,  B and  G'  into  solvability  condition,  which  is  given  by 


fO  ^^s,i)Coshk{h  + z) 
J-h  cosh  kh 


cosh  kh 


(A.29) 


in  which 


B{3,l)  = -iki4>^^’^)hx,= 


gh  1 

2nicoshfch 


for  slip  bottom  boundary  condition,  or 

for  non-slip  bottom  boundary  condition; 

/'°  ipf3,i)CoshA:(/i-|- 2)  h kiaix2  1 

J-h  cosh  kh  2cti  aj 

\ ^ hx2  . .ccg  . k^  - k^ 

and 


^f(3,l) 


^11  ,9  ig  9 

2tanhlfch^  I a I a + ar.K,  + gUax2  + 

+(|V/.-U-^U.V;,^i)a  + 0(e'‘) 


where 


cosh  4kh  -1-8  — 2 tanh^  kh 


8 sinh"*  kh 
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Substituting  those  expressions  into  solvability  condition  (A.29),  we  obtain  the  governing 
equation  combining  refraction  with  diffraction  and  current. 

2{ccgh  + <riU)a,  + 2a^Vay  - i{ccg  - V^)ayy  + 

+ iccg{k\  - k^)]a  + ikccgk'  | a ^ a - ■*  = 0 (A-30) 

cosh  kh  \ 2u) 

where  k — k —D.  The  last  term  in  the  equation  only  exists  for  non-slip  bottom  boundary 
condition.  It  is  exactly  the  same  as  we  have  developed  in  chapter  2. 


APPENDIX  B 

NON-SLIP  BOTTOM  BOUNDARY  CONDITION  APPLIED  IN  WAVE  EQUATION 


Appendix  B will  present  how  to  get  non-slip  boundary  condition  for  wave  equation  in 
Appendix  A. 

Following  Liu(1986),  we  assume  that  the  total  velocity  near  the  bottom  can  be  expressed 
as 


— tiyji  ”1“  ^ci  ^ (-2^  ”h  ^ S 


(B.l) 


where  u^,-  is  rotational  flow  to  satisfy  u,-  = 0 at  bottom,  subscripts  w and  c represent 
quantities  of  wave  and  current,  respectively,  and  S = is  boundary  layer  thickness. 

The  linearized  Navier-Stokes  equation  gives 

dui  1 dp  d^Ui 

~ dxjdxj 

If  there  is  no  rotational  flow,  (B.2)  becomes 

duti  1 dp  d^uti 

~ ~~p~d^i^^  dxjdxj 

where  uu  = u^i  + u^i 

Therefore  from  (B.2)  and  (B.3),  we  obtain 


dUri  _ ^ d^Uri 
dt  dxjdxj 

Assuming  u„-  is  periodic  in  time 


(B.4) 


results  in 


u„  = ui,e 


—iut 


— iojuii  = V 


dxjdxj 


(B.5) 
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where  uu  = uu{x,y,  z„),  in  which  z„  is  the  coordinate  perpendicular  to  the  bottom.  Defin- 
ing 


(B.5)  becomes 


1 


iuii  = 2 


d uii 

5^2 


(B.6) 


Solving  (B.6)  with  boundary  conditions 

z = 0 ^ cx) 


Uri  = 0 


z = -h  ^ = 0 u„-  = 

where  is  linear  wave  potential  as  defined  in  Appendix  A,  we  have 


d6 

1 r * 

1— ‘ “pI  2« 

-(,  + A)l 

(B.7) 

oy 

1 r (1  ~ * 

1.=--.  ‘M 

-{z  + h)] 

(B.8) 

From  the  continuity  equation,  we  have 
dWi 


dz 


^ dx  dy^ 

a2. 


^ r (1-0,  , l,d<t>dh  d<t>dh. 

+ ^)„-»exp|-^(,  + A)|  - 


26 

exp[  - (z  + ft)  ] 

After  integrating  with  respect  to  z , we  get 

^1  = —{i  + + ^^4)  •'Vhh  at  z = —h 

where  is  two  dimensional  gradient  operator. 

If  a.ssuming 

6 ~ 0(eO,  Vfc  ~ 0(eO 
Elquation  (B.IO)  at  order  of  O(e^)  becomes 


(B.9) 


(B.IO) 


(B.ll) 
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where  ^ heis  been  applied. 

Therefore,  we  get  finally,  at  bottom 


— + = 0 (B.12) 


which  is  at  order  of  O(e^). 


APPENDIX  C 

NUMERICAL  SCHEME  FOR  CONTINUITY-MOMENTUN  EQUATION 


Appendix  C will  give  the  general  numerical  model  scheme  for  continuity  equation  and 
momentum  equations. 

Following  the  linearized  derivation  in  chapter  2,  we  may  write  the  entire  equations  of 
(2.11),  (2.21)  and  (2.22)  in  chapter  2 as 


dri  dqx  ^<iv 

dt  dx  dy 

dqx  ^drj 

■ ''''  “ 


where 


qx  — U D and  —V  D are  total  flux  in  x and  y direction  respectively. 

= -^  \ u'f  \x  and  Cy  = ^ I |y  as  define  in  (2.136)  in  chapter  2. 
rrix  and  rriy  are  the  rest  term  in  x and  y momentum  equations,  respectively. 
Casting  (C.1-C.3)  in  difference  form,  we  obtain 

ff}+i  _ . j - 

’•'a,  - O + A^f’S^  - = ° 

+ m.,„-  = 0 

= 0 

where  the  quantities  without  time  levels  denote  the  values  at  time  level  n. 

In  the  matrix  form,  those  equations  become 

^ ^ + AtM  = 0 


(C.l) 

(C.2) 

(C.3) 


(C.4) 

(C.5) 

(C.6) 


(C.7) 


where  Sx  and  Sy  are  central  spatial  differential.  The  matrixes  W,  A,  B,  C and  M are  defined 
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( 

0 

1 

0 ^ 

( 0 

0 

1 A 

A = 

gD 

0 

0 

B = 0 

0 

0 

\ 

0 

0 

0 J 

V gD 

0 

oy 

(n  ] 

( ° 1 

( M 

w = 

Qx 

c = 

Cx 

M = 

rrix 

( J 

< J 

1 7 

where  is  transposed  matrix  of  C. 

Introducing  Px  — = CAt,  Eq.(C.7)  can  be  rewritten  as 

{l  + ^x  + ^y  + ^W-  AtM  (C.8) 

By  adding  -W)  and  ^cPy{W'^+^  -W)  into  above  equation,  both  are  high  order 

truncation  errors,  we  have 


(1  + + /?,)(!  + ^y)W^+^  = (1  + ^x^y  + ^,l3y)W  - AtM 


Eq.(C.9)  can  be  easily  split  into  two  parts 


(C.9) 


{l  + ^x  + ^c)W*  = {l-^y)W -AtM  (C.IO) 

+ =W*+^yW  (c.ll) 

(C.IO)  and  (C.ll)  emerge  into  (C.9)  once  W*  is  eliminated.  Writing  (C.IO)  and  (C.ll)  in 
detail,  results  in 


(1  + AtCx)q*x  + -^gDSxfj*  = qx-  Atrrix 
(1  + AtCy)q*  = qy  - ^gDSy^-  Atrriy 


(C.12) 

(C.13) 

(C.14) 


and 


= r + (C.15) 

= «:  (C.16) 

(C17) 
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Substituting  (C.16)  into  (C.12)  and  (C.13),  we  obtain  x-sweep 

I n*  ~n  + + ^Syqy  - 0 

[ (1  + -q^  + ^gDS^ri*  + Atm^  = 0 

and  after  combining  (C.14)  and  (C.17),  we  get  y-  sweep 

I r-^"-^*  + ^y^y<ir^-^ySy<ly  = 0 

[ (1  + + ^gD6yfr+'^  + Atmy  = 0 


(C.18) 


(C.19) 


In  the  second  equation  of  (C.19)  the  high  order  term  ^^^gDcySy{rf^+^  - ff)  has  been 
omitted. 

The  above  two  set  equations  are  almost  identical  to  those  developed  in  chapter  3. 


APPENDIX  D 

DOUBLE-SWEEP  METHOD  IN  SOLVING  NUMERICAL  MODELS 


Appendix  C will  present  the  details  of  double-sweep  method  used  in  chapter  3. 
Continuity-Momentum  Equations 
(1)  X-Sweep 

Elquation  (3.21)  in  chapter  3 read 


- Pi-ijurj-'  = c^j 

= Bij 


(D.l) 

(D.2) 


where 
Rl  = 
PiJ  = 


gAt 

Ax 

At 

2Ax 


(A+i,y  + Dij) 


At 


A.i  Vi,j  . [ -I-  Dij)  - Vij{Dij  + D,j_x)  ] 


— 1 -|- 
Bij  = Uij 


2Ay 
fAt 


{Dij  -f  Di^i  j) 

At 


M,j 


At 

At 


+^'  i,y)(A.y+i  Uij  i)  ('5rrx).-ij]  ^ 

[ (Ay)«,y+1  ~ (Ay)»J-l  + — (5j;y),_X,y_l  ]}  {F  Lyy)ij  + {F  Lyx)ij 

The  quantities  without  superscripts  denote  the  value  at  time  level  n.  | (ut),-,y  |x,  {FLyy)ij 
and  {FLyx)ij  are  defined  in  (3.18)  in  chapter  3. 

The  solutions  of  rff  j and  are  assumed  with  the  forms  as 

= Fli  + EliUpf^ 

= F2,.,  + 
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(D.3) 

(D.4) 
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Substituting  (D.3)  and  (D.4)  into  (D.2)  and  (D.l),  respectively,  we  obtian 


fl*ij  = Fli  + EliUpf^ 


(D.5) 


which  is  the  same  as  (3.23)  in  chapter  3.  The  definitions  of  FI,  El,  F2  and  E2  are  also 
given  in  chapter  3. 

(2)  Y-Sweep 

Equation  (3.22)  in  chapter  3 gave 


+ Qijvr/A  - = Zi.i  (D.6) 

= Sij  (D.7) 


where 


R2 

Qi,j 

ZiJ 

Fi.: 

SiJ 


g^t 

Ay 

At 


(A,i+i  + 
At 


/At 

^ + 2( + A-,y-i)*  ' 

+ Ui,j  + + C/.- ,_i)"+i(V,+i.,-  - V;_i,,) 

+ ('S'iv)t+i.;-i  - (5'iy),-ij_i]}  + -^{^vv)i,i  - ('S'vv)m-i]  + {FL^y)ij  + {FLyy)ij 


where  | (ut),j  |y,  (^FLxy)tj  S'Hd  (^FLyy)ij  are  defined  in  (3.20)  in  chapter  3.  The  quantities 
with  the  subscripts  * are  up-date  water  depth  defined  as  — h + ff* . 

As  done  in  x-sweep,  finally,  we  have 

[ = F3j  + ESjVp^J^^ 

1 (D.8) 


which  is  the  same  a.s  (3.24)  in  chapter  3. 
Wave  Equation 
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Rewriting  (3.28)  in  chapter  3 by  the  orders  of  grids,  we  obtain 


/\  j; 

A+i,y+if  — CO. 


4Ay“  ■ 2(Ay)2 


+ A>„|1  + 


4Ay 

Ai 

2(Ay) 


4Ay 


+i>i 


2^l«+i,y  ] — Qij 


(D.9) 


or 

-^^«+l.J-^+l,y+l  + ^^i+l,j  ■^+1,3  + = Qi,j  (D.IO) 

where 


B2  ■ = 

’■'  C2,.y  + C2._i,y 

53  ■ = g^M+i-<^3.,y-i 
C2,.y 

•^2,-  y = 1 + Mj  y 

Qij  is  the  all  terms  at  rows  below  to  t -f  1. 


~ ^^«>y^»,y+l  ^4,-,y-A,-,y  — F3ijAi  j-i 


in  which 


7^4, _y  — - 1 + M,_y 

After  some  manipulations  as  did  in  x-sweep  continuity- momentum  equation,  Eq.(D.lO) 
can  be  rewritten  as,  in  an  implicit  marching  approach 


A, — R2j  -(-  RlyA^- . 


(Dll) 


Rly  = 


-F3i 


.3+1 


F2.- y+i  + FI 

».y+i-^iy+i 


where 
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R2  ■ = 

^ + ■^It.j+i^ly+i 

The  double  sweep  will  start  from  offshore,  i=l,  in  y direction,  and  give  the  results  row  by 
row  till  shoreline,  i=idry. 
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